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This  study  focuses  on  state  space  partitioning  techniques  for  computing  mea-  ft 

sures  related  to  the  operation  of  stochastic  systems.  These  methods  iteratively  de¬ 
compose  the  system  state  space  until  the  measure  of  interest  has  been  determined. 

The  information  available  in  each  iteration  yields  lower  and  upper  bounds  on  this 
measure,  and  can  be  used  to  design  efficient  Monte  Carlo  estimation  routines.  We 
present  here  new  theoretical  results  identifying  strategies  for  significantly  enhancing 
the  performance  of  these  algorithms.  Using  these  results,  we  describe  a  generalized 
algorithm  that  can  easily  be  tailored  to  address  a  variety  of  problems.  We  next  use  ft 

this  algorithm  to  analyze  two  important  models  in  the  area  of  stochastic  network 
optimization.  The  first  concerns  the  probabilistic  behavior  of  minimum  spanning 
trees  in  graphs  with  discrete  random  arc  weights.  Sp<  :ifically,  we  compute  the  prob-  ^ 

ability  distribution  of  the  weight  of  a  minimum  spanning  tree  and  the  probability 
that  a  given  arc  is  on  a  minimum  spanning  tree.  Both  of  these  problems  are  shown 
to  be  #P-hard.  The  second  model  considers  minimum  cost  flows  in  networks  with 
discrete  random  arc  costs  and  capacities.  We  consider  the  case  of  statistically  inde-  ft 

pendent  costs  and  capacities  for  each  arc  as  well  as  the  case  in  which  the  cost  and 
capacity  of  each  arc  change  simultaneously.  In  each  case,  we  show  that  the  evalua¬ 
tion  of  the  distribution  of  the  minimum  cost  flow  for  a  fixed  demand  configuration 
is  a  #P-hard  problem.  Numerical  examples  are  given  throughout  the  thesis. 
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CHAPTER  1 


Introduction 


Much  work  has  been  done  on  the  various  problems  that  arise  in  the  area  of 
network  optimization,  resulting  in  a  rich  and  elegant  theory.  We  assume  familiarity 
with  this  theory,  and  will  use  its  terms  and  results  without  lengthy  explanation. 
See  [24], [28], [29]  for  an  overview  of  network  theory.  Network  optimization  problems 
are  unique  among  integer  and  combinatorial  optimization  problems  in  that  efficient 
algorithms  for  their  solution  generally  exist.  These  tools  can  solve  an  extensive  array 
of  optimization  problems,  including  many  that  do  not  naturally  involve  networks 
but  can  be  modeled  using  networks. 

In  practice,  unfortunately,  the  situation  is  complicated  by  the  fact  that  some  of 
the  network  parameters  (e.g.,  arc  weights  or  capacities)  either  are  not  be  known  with 
certainty  or  change  in  some  probabilistic  fashion.  In  such  cases,  it  would  likely  be 
necessary  to  make  multiple  optimization  runs  using  a  number  of  possible  values  for 
those  parameters.  In  fact,  in  order  to  best  gain  insight  into  an  uncertain  situation 
one  might  assign  probabilities  to  many  possible  parameter  settings  and  solve  each 
resulting  optimization  problem  separately.  The  resulting  distribution  would  then 
enable  one  to  make  probabilistic  statements  about  the  value  of  an  optimal  solution. 

Example  1,1  Suppose  a  communication  network  is  to  be  built  to  connect  the  nodes 
in  an  undirected  graph,  using  a  least-weight  subset  of  possible  arcs  ( weight  may 
represent  cost,  amount  of  material,  etc.).  During  the  design  phase,  the  weights 
of  these  arcs  are  not  known  with  certainty;  rather,  they  are  independent  discrete 


T*ble  1.1:  Input  Distribution  for  Example  1.1  and  Example  1.2. 
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Arc 

Weight  (Probability) 

1  =  0,2) 

2.0  (0.90) 

8.0  (0.08) 

12.0  (0.02) 

2  =  (1,3) 

10.0  (0.85) 

24.0  (0.12) 

35.0  (0.03) 

3  =  (1,4) 

6.0  (0.88) 

18.0  (0.10) 

24.0  (0.02) 

4  =  (2,5) 

17.0  (0.75) 

35.0  (0.20) 

50.0  (0.05) 

5  =  (2,6) 

3.0  (0.68) 

7.0  (0.25) 

10.0  (0.07) 

6  =  (2,3) 

12.0  (0.85) 

22.0  (0.11) 

30.0  (0.04) 

7  =  (3,7) 

10.0  (0.80) 

19.0  (0.14) 

24.0  (0.06) 

8  =  (3,8) 

4.0  (0.75) 

6.0  (0.17) 

10.0  (0.08) 

9  =  (3,4) 

3.0  (0.84) 

7.0  (0.13) 

12.0  (0.03) 

21.0  (0.85) 

33.0  (0.09) 

45.0  (0.06) 

11  =  (5,7) 

8.0  (0.77) 

14.0  (0.13) 

18.0  (0.10) 

15.0  (0.76) 

21.0  (0.22) 

25.0  (0.02) 

13  =  (3,6) 

5.0  (0.651 

10.0  (0.23) 

12.0  (0.12) 

14  =  (5,6) 

18.0  (0.94) 

27.0  (0.05) 

36.0  (0.01) 

15  =  (6  7) 

38.0  (0.12) 

50.0  (0.08) 

16  =  (7,8) 

30.0  (0.09) 

40.0  (0.04) 

17  =  (7,10) 

4.0  (0.75) 

9.0  (0.14) 

15.0  (0.11) 

(I 

00 
f  "  * 

9.0  (0.82) 

13.0  (0.15) 

20.0  (0.03) 

19  =  (8,9) 

15.0  (0.79) 

28.0  (0.15) 

45.. 0  (0.06) 

20  =  (7,9) 

10.0  (0.95) 

17.0  (0.03) 

30.0  (0.02) 

21  =  (9,10) 

8.0  (0.85) 

14.0  (0.10) 

25.0  (0.05) 

random  variables  with  the  probability  distributions  given  in  Table  1.1.  (For  example, 
arc  1  will  have  weight  2.0  with  probability  0.9,  8.0  with  probability  0.08,  etc.)  This 
is  a  stochastic  version  of  the  classic  minimum  spanning  tree  (MST)  problem.  If  we 
are  constrained  to  build  the  MST  at  a  weight  that  does  not  exceed  a  given  value  d, 
we  might  like  to  know  the  probability  that  we  will  be  successful.  In  fact,  it  would 
be  helpful  to  determine  the  entire  distribution  of  the  weight  of  a  MST  so  that  the 
expected  weight  and  higher  moments  may  be  computed. 

Example  1.2  Given  the  same  network  and  arc  weight  distribution,  we  are  inter¬ 
ested  in  computing  criticality  indices  of  arcs,  as  defined  below. 
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Definition  1.1  The  criticality  index  for  arc  e  is  the  probability  that  this  arc  belongs 
to  a  minimum  spanning  tree. 

The  arcs  with  the  largest  criticality  indices  are  most  likely  to  be  on  the  least-weight 
spanning  tree  described  above. 

Example  1.3  In  a  directed  flow  network  of  10  nodes  and  21  arcs  with  source  s  =  1, 
sink  t  =  10,  and  arcs  having  independent  random  costs  and  capacities  per  unit  of 
flow  transmitted  given  in  Tables  1.2  and  1.3,  respectively,  we  wish  to  derive  the 
distribution  of  the  cost  of  a  minimum  cost  s—t  flow  satisfying  a  demand  requirement 
for  15  units  at  node  10. 

Note  that  uncertainty  about  network  parameters  can  come  about  in  a  number 
of  ways.  In  the  examples  above,  we  had  the  type  of  uncertainty  that  is  necessarily 
involved  in  planning.  In  some  systems  the  random  variation  can  be  approximated 
by  discrete  distributions,  as  in  the  following  examples. 

Example  1.4  Data  packets  of  a  uniform  size  are  transmitted  over  the  links  of  a 
communication  network.  Each  packet  contains  information  represented  in  bits,  each 
of  which  is  subject  to  transmission  error  (with  probability  that  may  depend  on  the 
link  being  used).  From  this  we  derive  the  probability  that  a  packet  arrives  at  the 
other  end  of  a  given  link  containing  errors.  If  an  arriving  packet  contains  an  error, 
it  must  be  re-transmitted.  Thus  the  amount  of  time  it  takes  for  a  packet  to  be 
correctly  transmitted  across  a  link  depends  on  the  length  of  the  link  and  on  the  to¬ 
tal  number  of  transmissions  needed  for  successful  transmission  (a  geometric  random 
variable  if  we  assume  that  the  several  transmissions  of  a  packet  constitute  indepen¬ 
dent  Bernoulli  trials  with  constant  probability  of  success),  which  in  turn  depends 
on  the  particular  bit  error  characteristics  of  that  link.  Then  the  resulting  transmis¬ 
sion  time  distribution  is  naturally  discrete  and,  from  a  practical  standpoint,  can  be 
truncated  and  thus  be  considered  to  have  finite  state  space.  If  the  network  is  to  be 
connected  using  only  the  most  economical  links,  we  are  interested  in  investigating 
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Table  1.2:  Arc  Cost  Distributions  for  Example  1.3. 


Arc 

Arc  Cost  (Probability) 

1  =  (1.2) 

70  (0.5) 

73  (0.3) 

gFQHI 

2  =  (1,3) 

25  (0.5) 

35  (0.4) 

3  =  (1,4) 

42  (0.5) 

48  (0.3) 

fagtorol 

4  =  (2,5) 

26  (0.4) 

31  (0.3) 

55  (0.2) 

88  (0.1) 

5  =  (2,6) 

58  (0.4) 

70  (0.3) 

95  (0.3) 

6  =  (3,2) 

15  (0.6) 

73  (0.4) 

7  =  (3,7) 

65  (0.5) 

74  (0.4) 

75  (0.1) 

8  =  (3,8) 

59  (0.6) 

72  (0.3) 

98  (0.1) 

9  =  (4,3) 

21  (0.3) 

32  (0.3) 

85  (0.2) 

98  (0.2) 

89  (0.7) 

96  (0.3) 

11  =  (5,7) 

32  (0.6) 

48  (0.2) 

67  (0.2) 

63  (0.5) 

99  (0.5) 

13  =  (6,3) 

66  (0.8) 

85  (0.1) 

98  (0.1) 

14  =  (6,5) 

6  (0.4) 

15  (0.3) 

39  (0.2) 

58  (0.1) 

15  =  (6,7) 

2  (0.6) 

48  (0.4) 

16  =  (7,8) 

16  (0.3) 

18  (0.3) 

40  (0.2) 

52  (0.2) 

!G  (0.5) 

34  (0.4) 

71  (0.1) 

18  =  (8,4) 

96  (0.5) 

19  =  (8,9) 

17  (0.4) 

49  (0.4) 

53  (0.1) 

65  (0.1) 

20  =  (9,7, 

3  (0.5) 

30  (0.4) 

50  (0.1) 

21  =  (9,10) 

6  (0.5) 

12  (0.3) 

54  (0.1) 

66  (0.1) 

Table  1.3:  Arc  Capacity  Distributions  for  Example  1.3. 


Arc 

Arc  Capacity  (Probability) 

WM 

2  =  (1,3) 

BY  Dpj 

1MEI 

15  (0.6) 

3  =  (1,4) 

111 

BEBI 

4  =  (2,5) 

B  M 

10  (0.3) 

17  (0.5) 

5  =  (2,6) 

4  (0.1) 

7  (0.2) 

8  (0.3) 

10  (0.4) 

6  =  (3,2) 

7  (0.1) 

9  (0.2) 

15  (0.7) 

P 

CO 

II 

5  (0.1) 

8(0.1) 

10  (0.3) 

16  (0.5) 

8  =  (3,8) 

8(0.1) 

12  (0.9) 

9  =  (4,3) 

9  (0.1) 

13  (0.2) 

18  (0.7) 

10  =  (4,9) 

4  (0.2) 

7  (0.2) 

11  (0.2) 

113  (0.4) 

11  =  (5,7) 

8  (0.5) 

14  (0.5) 

12  =  (5,10) 

6  (0.1) 

15  (0.2) 

17  (0.7) 

13  =  (6,3) 

4(0.1) 

7  (0.2) 

12  (0.2) 

20  (0.5) 

14  =  (6,5) 

5  (0.3) 

12  (0.7) 

15  =  (6,7) 

5  (0.2) 

7  (0.3) 

8  (0.5) 

16  =  (7,8) 

5  (1.0) 

17  =  (7,10) 

6  (0.2) 

14  (0.8) 

18  =  (8,4) 

3  (0.1) 

5(0.1) 

9  (0.2) 

13  (0.6) 

19  =  (8,9) 

8  (0.2) 

10  (0.2) 

14  (0.6) 

20  =  (9,7) 

4  (0.1) 

7  (C  ?) 

9  (0.3) 

12  (0.4) 

21  =  (9,10) 

5  (0.3) 

12  (0.7) 

minimum  spanning  trees  using  the  probabilistic  structure  just  described.  (Local 
area  networks,  for  example,  often  use  a  tree  architecture.)  Criticality  indices  will 
identify  the  links  that  are  most  likely  to  be  economical. 

Example  1.5  A  given  number  of  trucks  are  available  to  transport  goods  from  city 
A  to  city  B  each  day.  Sometimes  (according  to  historical  probabilities  which  are 
known),  one  or  more  trucks  are  unavailable.  Thus  the  capacity  of  the  A-B  trans¬ 
portation  link  depends  on  the  number  of  available  trucks,  clearly  a  discrete  random 
variable  with  finite  state  space.  Investigating  minimum  cost  flows  in  a  network  made 
up  of  many  such  shipping  links  will  involve  discrete  random  capacities. 

Example  1.6  Often  an  excellent  model  for  randomly-occurring  phenomena  is  the 
Poisson  process.  Suppose  that  the  occurrence  of  some  delaying  phenomenon  along 
the  length  of  a  transportation  link  can  be  so  modeled.  Then  the  time  it  takes 
to  traverse  that  link  might  depend  on  the  number  of  such  phenomena  that  are 
encountered  while  traversing  the  link  and  the  time  for  traversing  that  link  is  clearly 
a  discrete  random  variable. 

The  cardinality  of  the  state  space  is  the  biggest  obstacle  in  the  evaluation  of 
probabilities  related  to  the  operation  of  networks  with  randomly  weighted  compo¬ 
nents.  For  instance,  the  network  in  Table  1.1  has  321  ~  10  billion  states  and  the 
evaluation  of  the  distribution  of  the  weight  of  a  MST  by  a  complete  enumeration  of 
the  state  space  requires  10  billion  executions  of  a  MST  algorithm.  If  each  arc  has 
four  possible  settings,  421  «  4  quadrillion  such  executions  are  needed.  The  state 
space  for  the  network  in  Example  1.3  has  nearly  1019  states.  It  is  this  exponential 
explosion  of  work  that  has  motivated  the  development  of  techniques  which  seek  to 
evaluate  only  a  fraction  of  these  possible  settings  without  sacrificing  the  quality  of 
the  resulting  information. 

In  1992,  Ball  et.  al.  surveyed  state  of  the  art  methods  for  evaluating  a  broad 
assortment  of  network  reliability  and  related  problems  [7].  In  that  paper  they  pro¬ 
posed  the  problems  in  Example  1.1  and  in  Example  1.3  as  multistate  performability 
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problems.  After  describing  the  stochastic  maximum  Sow,  shortest  path  and  PERT 
problems,  whose  solutions  they  treated  in  some  detail,  the  authors  pointed  out  that 
“more  elaborate  models  could  include  both  cost  and  capacity  information  (random 
or  deterministic)  and  include  a  wide  variety  of  measures  of  system  operation,”  with 
some  examples  being  the  stochastic  minimum  cost  Sow  and  minimum  spanning 
tree  problems.  They  found  no  satisfactory  solution  for  these  problems  and  sug¬ 
gested  classes  of  algorithms  to  be  applied  to  such  unsolved  problems.  They  stated 
that  exact  solution  methods  in  network  reliability  typically  involve  convolutions  or 
max(min)  operations  (arbitrary  series  of  which  are  known  to  be  #P-hard,  even 
for  variables  which  have  only  two  possible  values),  or  involve  transformations  into 
problems  with  two-state  variables  (a  process  that  typically  replaces  each  variable 
having  q  possible  values  by  q  two-state  variables).  Factoring,  in  which  a  state  space 
is  iteratively  partitioned  based  on  the  state  of  a  single  variable,  was  also  given  as 
an  important  solution  methodology.  A  variant  of  factoring  was  used  in  a  paper  by 
Doulliez  and  Jamoulle  [15]  that  solved  problems  related  to  maximum  flows  in  net¬ 
works  having  random  discrete  arc  capacities.  In  that  approach,  the  state  space  of 
arc  capacities  was  iteratively  decomposed  based  on  sets  of  state  points  into  known 
and  unknown  regions  until  the  probabilistic  behavior  of  all  points  in  the  state  space 
was  known.  The  survey  by  Ball  et.  al.  concluded  that  all  these  available  methods 
were  necessarily  inefficient  since  “essentially  all  reliability  problems  of  interest  are 
#P-complete.”  Thus  it  is  seen  that  our  proposed  problems  are  not  easy  ones. 

Some  special  cases  of  the  MST  problem  have  been  treated  in  some  detail.  Gilbert 
[18]  considered  the  problem  of  building  MSTs  to  connect  points  randomly  placed  in 
the  two-dimensional  unit  circle  according  to  a  Poisson  process.  Frieze  [17]  examined 
the  case  of  the  complete  graph  on  n  nodes  whose  arc  costs  are  given  as  independent, 
identically  distributed  random  variables.  Assuming  that  the  common  distribution 
function  satisfied  certain  smoothness  requirements,  he  derived  limiting  behavior  (as 
n  -♦  oo)  of  the  expected  MST  cost.  Steele  [31]  and  Bertsimas  [10]  strengthened  this 
result,  known  as  the  ((3)  limit.  Bertsimas  and  van  Rysin  [11]  developed  asymptotic 
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results  for  MSTs  built  to  connect  n  points  uniformly  and  independently  distributed 
in  the  (multi-dimensional)  unit  ball.  Kulkami  [23]  used  a  clever  Markov  process 
approach  to  solve  the  special  case  of  investigating  minimum  spanning  trees  when 
arc  weights  are  independent  exponential  random  variables.  Bailey  [5]  extended 
these  results  to  the  general  case  of  minimization  on  matroids  with  exponentially 
distributed  element  weights.  Jain  and  Mamer  [22]  studied  the  problem  of  finding 
the  mean  and  distribution  of  MST  cost  in  a  network  whose  arc  costs  are  independent 
random  variables.  They  developed  an  upper  bound  that  proved  to  be  better  than 
the  naive  bound  obtained  by  solving  the  deterministic  MST  with  expected  arc  costs. 
Procedures  for  approximating  the  distribution  function,  which  involved  convolutions 
and  Laplace-Stieltjes  transforms,  were  demonstrated  for  exponentially  distributed 
arc  costs. 

Thus,  satisfactory  solutions  to  the  MST  problems  of  Examples  1.1  and  1.2  do 
not  appear  in  the  literature.  Specifically,  no  results  deal  with  nonidentical  discrete 
arc  cost  random  variables  or  offer  procedures  applicable  here  for  exactly  determining 
the  distribution  of  the  MST  cost.  Furthermore,  results  relating  to  criticality  indices 
as  defined  here  in  a  MST  context  were  examined  only  in  [23]. 

Nor  is  the  problem  posed  in  Example  1 .3  addressed  in  its  full  generality  by  exist¬ 
ing  solutions.  Some  results  do  exist  for  a  few  special  cases,  as  we  now  discuss.  The 
literature  here  is  dominated  by  papers  dealing  with  various  forms  of  the  stochastic 
maximum  flow  problem.  The  case  in  which  only  arc  capacity  is  random  has  been 
treated  by  Doulliez  and  Jamoulle  [15],  Grimmett  and  Welsh  [19],  Fishman  and  Shaw 
[16],  Alexopoulos  and  Fishman  [3]  [4],  Alexopoulos  [1],  Nagamochi  and  Ibaraki  [27] 
and  others.  Chen  and  Zhou  [13]  considered  methods  for  system  enhancement  when 
only  demand  is  random.  Prekopa  and  Boros  [30]  evaluated  the  probability  that 
the  system  of  linear  inequalities  of  Gale  and  Hoffman  (which  provide  necessary  and 
sufficient  conditions  for  feasibility  of  demand)  are  satisfied.  This  work  assumed  ran¬ 
dom  capacities  and  demand.  Corea  and  Kulkami  [14]  used  first  passage  times  in  a 
particular  Markov  process  to  investigate  minimum  cost  transshipment  problems  in 
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networks  with  independent  exponential  arc  costs  and  no  capacity  limitations.  Fi¬ 
nally,  the  stochastic  shortest  path  problem,  where  axe  costs  are  random,  was  treated 
by  Alexopoulos  [2]. 

Thus  none  of  the  subject  problems  are  satisfactorily  treated  in  the  literature, 
yet  there  is  a  need  to  produce  the  kind  of  information  described  above.  Given 
the  prospects  of  dramatically  enlarging  the  state  space  with  a  transformation  or 
of  performing  exponential  numbers  of  convolutions  of  multistate  variables,  we  are 
motivated  by  the  factoring  approach  of  Doulliez  and  Jamoulle.  Similar  approaches 
have  been  used  to  solve  some  other  problems  involving  optimization  of  networks 
with  stochastic  components.  See  [3],  [4]  for  maximum  flow  problems  and  [2]  for 
shortest  path  problems.  Furthermore,  extensions  have  been  made  in  the  use  of 
information  obtained  during  a  state  space  decomposition  as  input  to  highly  efficient 
Monte  Carlo  sampling  routines  [1].  This  study  will  describe  a  generalized  version 
of  this  decomposition  approach,  and  then  use  its  expanded  capabilities  to  solve 
probabilistic  problems  related  to  minimum  spanning  trees  and  minimum  cost  flows. 
These  solutions  have  the  attractive  properties  that  exact  computations  typically 
require  small  numbers  of  iterations  relative  to  the  cardinality  of  the  state  space,  and 
that  the  algorithms  produce  upper  and  lower  bounds  for  the  measures  of  interest 
that  improve  after  each  iteration. 

Table  1.4  previews  a  few  results  from  Chapter  3  for  the  MST  model  in  Table  1.1 
(FORTRAN  77  programs  run  on  a  SUN  SPARCstation  IPC)  and  Table  1.5  has  rep¬ 
resentative  results  for  minimum  cost  flow  problems  (FORTRAN  77  code  executed 
on  a  SUN  SPARCstation  2).  Note  that  the  number  of  iterations  is  only  a  diminutive 
fraction  of  the  number  of  system  states.  The  minimum  spanning  tree  results  extend 
naturally  to  other  matroid  problems,  and  the  minimum  cost  flow  results  demon¬ 
strate  the  applicability  of  the  proposed  methods  to  a  variety  of  stochastic  network 
problems  (e.g.,  maximum  flow  problems,  shortest  path  problems,  assignment  prob¬ 
lems,  matching  problems,  etc.)  that  can  be  reduced  to  minimum  cost  flow  problems. 
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Table  1.4:  Sample  of  decomposition  results  for  Examples  1.1  and  1.2. 


Information  computed  exactly 

#  iterations 

CPU  time  (sec)  | 

P{MST  cost  <  60}  =  0.8267496 

1421 

P{MST  cost  <  90}  =  0.9999951 

127,683 

Entire  distribution  of  MST  cost 

634,405 

445.91 

P{arc  1  is  on  a  MST}  =  0.9393247 

12 

0.71 

P{arc  5  is  on  a  MST}  =  0.8227428 

200 

1.18 

P{arc  14  is  on  a  MST}  =  0.0099071 

169 

1.00 

Table  1.5:  Sample  of  decomposition  results  for  Example  1.3. 


Information  computed  exactly 

#  iterations 

CPU  time  (sec) 

P{MCF  cost  <  1900}  =  0.30378595 
P{MCF  cost  <  3500}  =  0.99999962 

78,488 

132,944 

310.25 

443.50 

In  Chapter  2  we  begin  by  describing  the  General  State  Space  Decomposition 
algorithm  and  point  out  its  utility  in  conducting  probabilistic  analysis,  not  nec¬ 
essarily  in  the  area  of  network  optimization.  In  Chapter  3  we  propose  a  number 
of  different  algorithms  for  stochastic  minimum  spanning  tree  problems  and  discuss 
extensions  to  general  matroid  settings.  Chapter  4  contains  algorithms  for  stochas¬ 
tic  minimum  cost  flow  problems.  We  conclude  in  Chapter  5  with  topics  for  future 
research.  Numerical  examples  are  given  throughout  this  study. 


CHAPTER  2 


General  State  Space 
Decomposition  Algorithm 


2 . 1  Int  roduct  ion 

In  this  chapter  we  describe  methods  for  evaluating  probabilities  related  to  sys¬ 
tems  of  independent  discrete  random  variables.  Faced  with  such  a  task,  one  must 
devise  ways  to  combat  the  associated  computational  intractability  mentioned  in 
Chapter  1.  Solution  methods  involving  state  enumeration  are  clearly  out  of  the 
question  since  the  state  spaces  of  these  systems  grow  exponentially  with  the  num¬ 
ber  of  variables. 

One  possible  approach  is  Monte  Carlo  simulation.  In  crude  Monte  Carlo,  state 
points  are  sampled  according  to  the  probability  distributions  of  the  input  random 
variables.  The  proportion  of  sampled  points  having  a  desired  property  is  an  estimate 
of  the  probability  that  the  modeled  system  has  that  property.  Unfortunately,  in 
order  to  have  an  estimate  with  acceptably  small  variance,  an  astronomical  number 
of  samples  must  be  drawn. 

Much  work  in  reliability  analysis  (a  special  case  of  the  present  setting  in  which 
system  components  are  given  by  binary  absence/presence  variables)  has  used  a  fac¬ 
toring  (or,  pivotal  decomposition)  approach,  in  which  the  state  space  is  iteratively 
conditioned  and  partitioned  based  on  the  state  of  a  single  component  [8).  Basically, 
factoring  iteratively  makes  use  of  the  following  fact  (where  the  event  0  corresponds 


to  the  system  operating  and  the  event  0,-  corresponds  to  component  i  operating): 

p{o)  *  pwmpm+p&mpm. 

These  procedures  have  had  good  success. 

When  the  random  variables  describing  the  system  are  allowed  to  take  on  arbi¬ 
trary  values,  though,  performance  can  suffer.  This  situation  can  grow  still  worse  if 
we  allow  the  variables  to  have  arbitrary  (finite)  numbers  of  possible  values  (as  we 
shall  do  in  the  present  general  setting).  Though  factoring  in  this  setting  has  been 
used  effectively  in  some  cases  ([21],  for  example),  in  general  it  is  only  effective  if  the 
system  being  studied  has  a  special  structure  that  can  be  exploited. 

In  1972,  in  a  paper  dealing  with  problems  related  to  stochastic  maximum  flows, 
Doulliez  and  Jamoulle  [15]  used  an  adaptation  of  the  factoring  approach  to  classify 
and  decompose  the  state  space  based  on  sets  of  state  points.  (The  connection 
between  factoring  and  our  algorithm  will  be  explained  in  Section  2.5.)  Since  then, 
similar  approaches  have  been  used  in  other  stochastic  network  settings  [1],[2],[3],[4]. 
The  purpose  of  this  chapter  is  to  describe  what  we  shall  refer  to  as  the  General  State 
Space  Decomposition  (GSD)  algorithm,  with  the  following  goals  in  mind: 

1)  To  define  terms  necessary  to  build  a  general  framework  within  which  we  may 
discuss  these  methods; 

2)  To  describe  a  general  input  problem  broad  enough  to  accomodate  all  known 
applications  of  the  Doulliez-Jamoulle  decomposition  approach  as  well  as  new 
classes  of  problems  both  within  and  without  the  realm  of  network  analysis; 

3)  To  describe  GSD,  which  should  not  only  expand  existing  methods  to  include 
the  broader  problem  class,  but  should  also  enhance  them  in  ways  that  can 
lead  to  more  efficiently  computed  bounds  on  the  distributions  of  interest;  and 

4)  To  derive  theoretical  results  comparing  GSD  to  alternative  decomposition 
strategies. 


We  shall  see  that  GSD  has  a  number  of  important  advantages  over  alternative 
methods.  (1)  Upper  and  lower  bounds  on  the  probability  of  interest  are  available 
after  each  iteration,  and  these  bounds  are  progressively  tighter.  (2)  The  information 
available  in  each  iteration  can  be  used  for  constructing  efficient  Monte  Carlo  pro¬ 
cedures  based  on  importance  and  stratified  sampling.  (3)  Sensitivity  analysis  can 
be  accomplished  very  easily  when  input  probabilities  change,  with  no  additional 
decomposition  effort.  (4)  GSD  is  maximally  efficient  (with  respect  to  minimizing 
the  number  of  partitions  formed  per  iteration)  among  a  large  class  of  methods. 

We  begin  in  Section  2.2  with  some  definitions  and  mathematical  preliminaries 
leading  up  to  a  general  problem  statement  in  Section  2.3.  Section  2.4  contains  the 
GSD  algorithm.  Section  2.5  establishes  theoretical  results  concerning  the  partition¬ 
ing  of  the  state  space.  Section  2.6  discusses  the  use  of  GSD  to  produce  bounds,  and 
Section  2.7  shows  how  these  bounds  can  be  used  as  inputs  to  Monte  Carlo  estimation 
routines.  Concluding  remarks  are  given  in  Section  2.8. 

2.2  Definitions  and  Preliminaries 

Consider  a  vector  X  =  (X\ , . . . ,  Xa)  of  discrete  random  variables  with  finite  state 
spaces.  Each  component  X}  takes  on  finite  values  u);(l)  <  ■••  <  u>j(rij)  (which  we 
refer  to  as  weights)  with  respective  probabilities  pj(l),. . . , pj{n}).  An  element  of  fl, 
the  state  space  of  X,  is  of  the  form  x  =  (tuj(/j),. . .  ,wa(la)),  where  lj  6  {1, . . .  ,rij} 
for  all  j.  In  the  interest  of  readability,  we  shall  often  refer  to  the  state  point  x  as 
simply  l  =  (I*,. . . , /„),  letting  the  index  lj  stand  for  Wj(lj).  Let  A  =  {1, . . .  ,a). 

Unless  otherwise  stated,  we  assume  that  the  random  variables  X}  are  mutually 
independent,  so  that  the  probability  of  the  state  point  x  is  calculated  easily  as 

P(x)  =  P{Xi=wl(ll),...,Xa  =  wa(la)}  (2.1) 

=  fl  ?>(*;)• 

;= i 

Often  we  will  be  dealing  not  with  individual  state  points,  but  with  special  types 
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Figure  2.1:  Discrete  interval  given  in  Example  2.1. 

of  sets  which  we  call  discrete  intervals. 

Definition  2.1  A  subset  of  fi  is  a  discrete  ( multi- dimensional)  interval  with  end¬ 
points  a  and  0  (a  <  0)  if  it  consists  precisely  of  all  state  points  l  —  (/j,...,/a) 
satisfying  /_,  €  {oj, . . .  ,0,}  for  all  j.  We  denote  this  interval  by  [a,/?]. 

Example  2.1  Let  a  =  2,  =  4,  and  n2  =  3.  Then  there  are  4  x  3  =  12  points 

in  fl,  as  shown  in  Figure  2.1.  The  state  points  shown  as  open  circles  constitute  the 
interval  [a,  0]  with  endpoints  a  =  (2,  l)  and  0  =  (3,3). 

The  probability  of  an  interval  can  be  easily  calculated,  again  due  to  independence 
of  the  random  variables,  as 


» 


> 


» 


i- 1  '~ai  | 

This  ease  of  calculation,  and  the  fact  that  intervals  can  be  stored  using  just 
the  endpoints  a  and  0,  make  intervals  very  attractive.  GSD  deals  exclusively  with 
intervals. 

Suppose  that  the  vector  X  represents  the  states  of  the  components  of  a  system 
which  operates  when  a  set  of  structural  constraints  S(X)  is  satisfied.  Vfe  further 
assume  that  the  validity  of  these  constraints  (i.e.,  whether  or  not  these  constraints  I 

are  satisfied)  can  be  checked  for  any  realization  of  X. 


14 


I 


The  following  definitions  set  the  stage  for  the  formal  description  of  GSD. 

Definition  2.2  Any  point  l  in  the  state  space  0  which  satisfies  all  constraints  in 

the  set  S(-)  is  a  feasible  point.  If  /  does  not  satisfy  S(-),  then  /  is  an  infeasible  point. 

< 

Definition  2.3  A  set  whose  states  have  been  classified  as  feasible  (infeasible)  is 
called  feasible  ( infeasible ).  Sets  with  unclassified  states  will  be  called  undetermined. 

Definition  2.4  The  maximal  feasible  (infeasible)  subset  of  Cl  is  called  the  feasible 
( infeasible )  region  and  is  denoted  ftF  (ft7).  Clearly,  ft  is  partitioned  as  ft  =  ftFUft7. 

The  purpose  of  our  algorithm  will  be  to  compute  the  probability  th&t  the  system 

y 

operates  (i.e.,  that  the  system  is  in  a  feasible  state),  denoted  PF  =  P{ftF}.  We 
shall  refer  to  this  measure  as  the  system  reliability. 

Thu3  far,  we  have  not  restricted  in  any  way  the  relationship  between  5(X)  and 
the  individual  random  variables  Xj.  We  must  now  do  so  by  requiring  a  sort  of 
monotonicity  of  the  satisfaction  of  S(X)  with  respect  to  each  Xj.  Formally,  we 
require  that  it  be  possible  to  unambiguously  classify  each  Xj  as  either  upper  feasible 
or  lower  feasible,  as  defined  below. 

Definition  2.5  X3  will  be  called  lower  feasible  if  for  every  realization  of  Xi,  i  f  j, 
all  levels  of  Xj  corresponding  to  feasible  state  points  are  less  than  all  levels  of  Xj 
corresponding  to  infeasible  state  points. 

* 

> 

Definition  2.6  Xj  will  be  called  upper  feasible  if  for  all  realizations  of  Xi,  i  f  j, 
all  levels  of  Xj  corresponding  to  feasible  state  points  are  greater  than  all  levels  of 
Xj  corresponding  to  infeasible  state  points. 

This  classification  induces  a  partition  of  the  components  of  X.  Henceforth,  we 
assume  that  the  random  variables  are  numbered  so  that  lower  feasible  variables 
are  numbered  as  X^, . . . ,  Xm  while  the  upper  feasible  variables  are  numbered  as 
Xm+\,  -  •  • ,  Xa.  We  shall  often  emphasize  this  distinction  by  writing  a  state  /  in 
partitioned  form  as 


^  (^1 » •  •  •  j  |^m+l  i  ■  •  ■  >  ^a)- 


Imposing  the  above  restriction  on  X  and  S(X)  was  necessary  to  establish  the 
following  useful  result. 

Proposition  2.1  Let  X  be  a  vector  of  a  indepedent  discrete  random  variables  with 
variables  X\  through  Xm  being  lower  feasible  and  variables  Xm+i  through  Xa  being 
upper  feasible.  Then  for  all  l  £  Cl: 

(i)  If  l  is  feasible,  then  the  set 

F  —  {l  £  Cl  :  1  <  U  <li  for  i  —  1, . . . ,  m 

h  <  h  <  n«  for  i  —  m  +  1 , . . . ,  a} 

is  ftosible. 

(ii)  If  l  is  infeasible,  then  the  set 

I  =  {l  £  Cl  :  /,  <  h  <  n,  for  i  =  1 , . . . ,  m 
1  <  /<  <  /i  for  i  =  m  +  1 , . . . ,  a} 

is  infeasible. 

Proposition  2.1  is  valuable  because  it  allows  us  to  classify  entire  intervals  in  fl  as 
feasible  or  infeasible  based  solely  on  the  feasibility  of  a  single  point.  It  also  polarizes 
each  interval  so  that  it  has  a  “feasible  end”  and  an  “infeasible  end.” 

Definition  2.7  Consider  an  interval  [a,  0)  in  0.  The  point 

iQ\0)  (®1)  ■  •  ■  t  |/^m  + 1 1  i  0a  ) 

will  be  called  the  extreme  feasible  candidati  ;  ‘hr  point 

(0\o)  =  (0\y - 0m1'  ■  .1.) 
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will  be  called  the  extreme  infeasible  candidate. 


Note  that  this  definition  does  not  imply  that  the  extreme  feasible  candidate  is  a 
feasible  point  or  that  the  extreme  infeasible  candidate  is  an  infeasible  point.  If  the 
extreme  feasible  candidate  is  found  to  be  feasible,  we  may  refer  to  it  as  the  feasible 
extreme.  Likewise  for  the  infeasible  extreme.  Indeed,  Proposition  2.1  implies  that  if 
is  infeasible  then  all  of  [a,  /3\  is  infeasible.  Similarly,  if  (0\ a)  is  feasible  then 
so  is  all  of  [a,  f$\.  Thus  Definition  2.7,  in  conjunction  with  Proposition  2.1,  hints 
at  a  key  idea  of  GSD:  First  try  to  classify  an  interval  as  feasible  or  infeasible  by 
evaluating  its  extremes.  If  this  fails,  start  at  the  feasible  (or  infeasible)  extreme 
and  push  in  as  far  as  possible  while  remaining  feasible  (infeasible),  yielding  a  large 
feasible  (infeasible)  interval.  Following  the  general  problem  statement  in  Section  2.3, 
we  shall  describe  in  Section  2.4  the  manner  in  which  this  idea  is  put  to  use. 

Before  getting  into  the  technical  details  of  GSD,  a  comment  concerning  system 
modeling  is  in  order.  Without  significant  modification  GSD  has  the  capability  for 
limited  representation  of  simple  dependencies  among  variables.  Specifically,  some 
of  the  variables  Xt, . . . ,  Xa  can  represent  vectors  of  dependent  variables.  This  tech¬ 
nique  can  result  in  more  realistic  modeling  (as  opposed  to  treating  these  attributes 
as  independent)  and  at  the  same  time  reduces  the  cardinality  of  the  state  space. 
(For  example,  replacing  3  variables,  each  with  4  levels,  by  a  triple  having  4  levels 
reduces  the  cardinality  of  fl  by  a  factor  of  42  =  16.)  Of  course,  the  realizations  of  Xj 
must  be  chosen  and  ordered  so  that  Xj  can  be  classified  as  upper  or  lower  feasible. 
For  example,  in  Chapter  4  we  shall  use  Xj  =  ( cost  of  arc  j,  capacity  oj  arc  j)  in 
the  stochastic  minimum  cost  flow  setting.  By  themselves,  cost  is  lower  feasible  and 
capacity  is  upper  feasible.  As  an  ordered  pair,  then,  high  capacity  must  be  paired 
with  low  cost  so  that  they  may  be  jointly  classified.  Fortunately,  it  is  reasonable 
in  many  applications  that  units  of  scarce  capacity  are  expensive  or  cause  extended 
delays. 


2.3  General  Problem  Statement 


» 


Recapitulating,  we  have  a  stochastic  system  represented  by  a  vector  X  = 
{Xu...  ,Xm\Xm+i, . . .  ,Xa)  of  a  mutually  independent  discrete  random  variables 
with  finite  support.  We  wish  to  calculate  the  system  reliability,  defined  as  the  prob¬ 
ability  that  the  constraints  S(-)  are  satisfied  (i.e.,  the  probability  of  the  feasible 
region  ftF).  We  assume  that  the  variables  X\  through  Xm  are  lower  feasible  while 
the  variables  Xm+i  through  Xa  are  upper  feasible.  This  causes  ft  to  be  ordered 
in  the  sense  of  Proposition  2.1.  We  proceed  to  describe  the  General  State  Space 
Decomposition  algorithm. 

2.4  The  GSD  Algorithm 

We  are  now  prepared  to  discuss  a  typical  iteration  of  GSD,  in  which  an  undeter¬ 
mined  set  is  partitioned  into  (non-overlapping)  feasible,  infeasible  and  undetermined 
intervals.  The  probabilities  of  the  feasible  intervals  are  accumulated  as  a  lower  bound 
for  PF,  the  infeasible  intervals  are  discarded,  and  the  undetermined  intervals  are 
filed  in  a  list  as  input  to  subsequent  iterations.  Following  the  description,  in  Sec¬ 
tion  2.5  we  present  theoretical  results  showing  the  performance  of  this  algorithm  to 
be  optimal  in  terms  of  producing  fewest  new  undetermined  intervals  per  iteration. 

There  are  two  basic  approaches.  Decomposition  from  the  feasible  extreme  point 
(or,  feasible  decomposition)  begins  at  the  feasible  extreme  and  attempts  to  push  in 
as  far  as  possible,  producing  one  large  feasible  interval.  Decomposition  from  the 
extreme  infeasible  point  (or,  infeasible  decomposition),  on  the  other  hand,  drives  in 
from  the  infeasible  extreme  to  remove  one  large  infeasible  interval.  We  shall  describe 
feasible  decomposition  in  detail:  infeasible  decomposition  proceeds  symmetrically. 

We  begin  with  an  undetermined  set  U  =  [a,  0],  As  a  preliminary  screening, 
evaluations  are  made  at  both  extremes  of  the  set.  If  S(-)  is  not  satisfied  at  (a|/?)  (that 
is,  if  the  extreme  feasible  candidate  is  an  infeasible  point),  then  by  Proposition  2.1 
the  entire  set  is  infeasible,  it  is  discarded  and  the  iteration  is  ended.  Otherwise, 
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the  state  (a|/?)  is  feasible.  If  (/?|a)  also  satisfies  S(-),  then  by  Proposition  2.1  U 
is  feasible,  the  probability  of  U  is  added  to  the  accumulated  lower  bound  for  PF, 
and  again  the  iteration  is  ended.  If  not,  then  U  clearly  contains  both  feasible  and 
infeasible  points  and  we  proceed. 

The  major  step  of  feasible  decomposition  is  the  identification  of  a  feasible  set. 
The  exact  mechanics  of  this  step  depend  on  the  particular  problem  being  solved. 
We  seek  a  feasible  point  l  which  is  easily  derived  and  which  will  yield  as  large  a 
feasible  set  as  possible.  We  begin  with  our  only  known  feasible  point,  (a|/3).  We 
then  increase  the  weights  of  the  lower  feasible  components  and  decrease  those  of  the 
upper  feasible  components  as  much  as  possible  while  remaining  feasible. 

In  practice,  we  generally  stop  short  of  “as  much  as  possible”  in  the  interest  of 
efficiency.  Often  a  fairly  deep  cut-off  state  can  be  derived  with  reasonable  effort, 
and  to  push  a  moderate  amount  further  would  incur  significant  cost.  Additionally, 
in  many  cases  it  may  be  difficult  even  to  determine  the  extent  to  which  further 
pushing  is  possible.  Fortunately  this  is  not  necessary,  for  Proposition  2.1  guarantees 
that  any  feasible  point  yields  a  feasible  interval.  Just  how  much  can  be  removed  at 
a  time  will  likely  depend  on  the  problem  being  solved. 

After  completing  this  step,  we  obtain  a  feasible  interval 


F  =  {/  6  U  :  otj  <  lj  <  lj  for  j  =  1, . . . , m 

h  <  lj  <  Pj  for  j  =  m  +  1, ...  ,a) 


(2.3) 


and  add  the  probability  of  F 


P{F)  = 


m  h 

n  e  Pj(i) 

j=l  i-a3 


a  Pj 

n  i >j(*) 

[j=m+l  ,=lj 


(2.4) 


to  the  lower  bound  for  PF. 
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We  can  stop  after  determining  F,  breaking  up  the  remainder  of  U  into  up  to 
a  undetermined  intervals  as  we  indicate  below.  In  some  applications  of  GSD,  we 
shall  do  precisely  this.  However,  since  we  wish  to  account  for  as  much  undetermined 
probability  as  possible  in  each  iteration,  it  makes  sense  to  try  to  remove  infeasible 
sets  as  well.  This  cannot  be  done  naively.  As  we  shall  soon  see,  blindly  removing 
even  one  infeasible  set  in  addition  to  F  can  increase  the  number  of  new  undetermined 
sets  by  as  much  as  a.  The  following  technique  removes  up  to  a  infeasible  sets  in 
addition  to  F  and  cannot  increase  the  number  of  undetermined  sets  already  needed 
to  partition  U\F.  (In  fact,  the  number  can  be  reduced  since  some  of  the  sets  in  the 
partition  of  U\F  can  be  found  to  be  infeasible.) 

In  order  to  identify  infeasible  sets,  we  turn  again  to  Proposition  2.1.  This  step 
is  carried  out  individually  and  independently  for  each  of  the  a  components.  For 
component  j,  this  step  consists  of  starting  at  the  extreme  feasible  point  (a|/3)  and 
pushing  only  this  component  (increasing  its  weight  if  it  is  lower  feasible,  decreasing 
its  weight  if  it  is  upper  feasible)  as  far  as  possible  while  remaining  feasible.  Here  it 
is  essential  that  we  push  as  far  as  possible. 

Definition  2.8  We  define  the  limiting  feasible  index  for  component  j  by 
lj  =  max{7  :  <  7  <  /?,,  (oti, . . .  ,aj_i,7,orJ+1, . . .  ,am\0)  is  feasible} 

if  Xj  is  lower  feasible,  and  by 

lj  =  min{7  :  aj  <  7  <  #,  (ar|0TO+1, . . .  7, ft+i,  •  •  •  ,A»)  is  feasible} 

if  Xj  is  upper  feasible. 

Note  that  for  lower  feasible  variables  j  we  have  lj  >  lj  so  that  if  lj  =  (i}  we  can  set 
lj  =  0j  with  no  additional  effort.  Similarly,  for  upper  feasible  variables  j  we  have 
lj  <  lj,  allowing  us  to  set  /,  =  a,  if  lj  =  ar 

If  we  are  able  to  push  lj  all  the  way  (to  /?,  if  X:  is  lower  feasible,  to  a3  if  X3 
is  upper  feasible),  no  infeasible  set  is  produced.  Otherwise,  we  know  that  pushing 
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the  weight  of  component  j  one  more  level  while  leaving  all  other  variables  at  (a\P)  ' 

yields  an  infeasible  point,  and  hence  an  infeasible  set.  If  Xj  is  lower  feasible  (that 
is,  if  j  <  m),  the  infeasible  set  is 

Ij  =  {leU:  4  +  1  <lj<Pj  » 

<*i  <  li  <  Pi  for  t  ^  j},  (2.5) 

while  if  Xj  is  upper  feasible  (j  >  m  +  1)  this  set  is 

/j  =  {/  €  U  :  <*j  <  lj  <  lj  -  1  1 

«.  <  h  <  Pi  for  t  /  j).  (2.6) 

Note  that  the  sets  lj  are  not  disjoint  (indeed,  aU  contain  the  point  (/?)«)).  In  most 
cases  this  is  not  a  problem  since  the  infeasible  sets  are  discarded  en  masse.  In  some 
applications,  though,  infeasible  sets  are  retained  for  further  analysis,  as  we  discuss 
later.  In  those  cases  it  is  necessary  to  partition  the  infeasible  sets.  Recall  that  the 
union  of  non-disjoint  sets  Au . . . ,  Aa  can  be  partitioned  into  the  (disjoint)  sets  ) 

Ax  =  Ax 
A2  — 

^3  =  U  A2)  & 


Aa  =  Aa\(Ax  u---uia_i). 

Applying  this  principle  to  the  present  setting,  we  obtain  new  infeasible  sets  lj  as 
follows.  If  Xj  is  lower  feasible,  then 

lj  =  { /  €  U  :  a,  <  /,  <  li  for  i  <  j 

b  +  1  <  h  <  Pi  (2.7) 

«.  <  U  <  Pi  for  i  >  j}, 

whereas  if  X}  is  upper  feasible,  then 

lj  =  { /  €  U  :  a<  <  li  <  li  for  i  <  m 


» 


I 


» 
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I 


li  <li<  fa  for  m  <  i  <  j 
a,  <lj<ij-  1 
o,  <  /,•  <  Pi  for  t  >  j). 


(2.8) 


ft 


ft 


Notice  that  the  resulting  infeasible  sets  are  intervals. 

We  have  thus  found  one  feasible  set  (by  finding  the  feasible  point  l)  and  up 
to  a  infeasible  sets  (by  finding  the  limiting  feasible  index  lj  for  each  j).  All  that 
remains  is  to  account  for  the  portion  of  U  that  remains  undetermined.  We  wish 
to  represent  U\(F  U  (U“=1/j))  as  a  disjoint  union  of  intervals.  We  account  for  the 
points  “between”  F  and  lj  for  each  j  by  defining  the  set 


Uj  —  {l  £  U  :  lj  +  1  <  l,  <  lj 
Of.  <  U  <  h 
h  <  ii  <  Pi 


for  i  j,i  <  m 

for  i  ^  j,  i  >  m}  ,  if  j  <  m 


=  {l€U:  <  />  <  fi  -  1 


A 

cti  <  li  <  li  for  i  ^  j,  i  <  m 

h  <  li  S  Pi  for  i  /  j,  i  >  m}  ,  if  j  >  m. 


(2.9) 


To  verify  this,  note  that  putting  component  j  beyond  l}  is  sufficient  to  yield  points 
outside  of  F  independent  of  the  levels  of  the  other  components.  Also,  no  component 
is  permittea  .o  be  at  a  level  beyond  the  index  /  since  this  would  yield  points  already 
accounted  for  by  U“=1/j. 

Unfortunately,  the  undetermined  sets  Uj  overlap  and  their  union  must  be  parti¬ 
tioned.  This  is  done  by  defining  the  intervals 


ft 


ft 


ft 


ft 
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ft 
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Uj  =  {/  €  U  :  a,-  <  U  <  U  for  i  <  j 

lj  +  l<lj<  h 

ft 

or,  <  U  <h  for  j  <  i  <  m 
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U  <  h  <  A*  for  *  >  to)  ,  if  j  <  m 


(2.10) 


a=  {/  £U  :  a,  <  li  <li  for  i  <  m 

5  <  A  <  A  for  m  <  *  <  j 

A  <  A  <  A  for  *  >  i)  ,  if  j  >  m. 

It  is  now  easy  to  show  that  U  =  FU  {U“=1/,}  U  {U*_lI/,},  thus  completing  the 
partition  of  (7.  If  we  chose  not  to  remove  infeasible  sets,  so  that  limiting  indices 
were  not  found,  we  replace  /,  in  this  description  by  A  if  i  <  m,  or  by  a,  if  i  >  m.  We 
calculate  P{F)  using  (2.4)  and  add  the  result  to  PtF,  the  accumulated  lower  bound 
for  PF.  The  undetermined  sets  are  filed  in  their  list,  and  the  infeasible  sets  are 
discarded  (or  retained  for  further  analysis).  The  iteration  is  over.  An  undetermined 
interval  is  then  selected  from  the  list  and  the  next  iteration  begins. 

To  summarize,  procedure  FEASIBLE(fZ)  describes  the  decomposition  of  a  sin¬ 
gle  undetermined  set  U  =  (a,  /?].  The  list  U  contains  all  remaining  undetermined 
intervals. 


Procedure  FEASIBLE(C7) 

Step  1  Start  with  the  state  (a||d)  and,  while  S(-)  remains  satisfied,  determine  a 
feasible  state  l  by  pushing  up  from  a j  for  j  <  m  and  down  from  A  for  j  >  m. 

Step  2  Determine  limiting  feasible  indices  lj  as  in  Definition  2.8.  (If  infeasible  sets 
are  not  to  be  removed,  set  /  =  (/9|a).) 

Step  3  Set  Pf  =  PF  -f  P{F},  where  P{F }  is  calculated  as  in  (2.4). 

Step  4  For  j  =  1, . . .  ,a:  If  lj  ^  lj,  file  Uj  in  If,  where  Uj  is  defined  by  (2.10). 

As  indicated  previously,  decomposition  of  a  given  interval  can  proceed  from  either 
the  extreme  feasible  point  or  the  extreme  infeasible  point.  Infeasible  decomposition 
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is  completely  symmnetric  to  feasible  decomposition.  We  find  an  infeasible  point  l 
by  starting  from  (/?|a),  and  define  the  infeasible  interval 

I  =  {leU:  lj<lj<Pj  for  j  =  1 , . . . ,  m  (2.11) 

otj  <lj<lj  for  j  =  m  +  1 , . . . ,  o}. 

If  we  also  desire  to  remove  feasible  sets,  we  continue  as  follows. 

Definition  2.9  We  define  the  limiting  infeasible  index  for  component  j  by 

lj  =  min{7  :  otj  <  7  <  fa,  (Pi,...,Pj„l,'i,pj+l,...,Pm\a)  is  infeasible}, 

if  Xj  is  lower  feasible,  and  by 

lj  =  max{7  :  otj  <  7  <  (p\am+i, . . .  ,«j_i,7,aJ+i,.. .  ,a0)  is  infeasible} 
if  Xj  is  upper  feasible. 

Note  that  in  infeasible  decomposition  we  have  l:  <  lj  for  j  <  m  and  lj  >  lj  for 
j  >  m. 

The  indices  lj  define  up  to  a  disjoint  feasible  intervals  F}.  If  X}  is  lower  feasible 
(i.e.,  if  j  <  m),  then 

Fj  =  {/  €  V  :  li  <  li  <  0i  for  i  <  j 

otj  <  lj  <lj-l  (2.12) 

0.  <  h  <  Pi  for  i  >  j}, 

whereas  if  Xj  is  upper  feasible  (j  >  m),  then 

Fj  =  {/  €  U  :  U  <  li  <  Pi  for  i  <  m 

a,  <  /,•  <  /,  for  m  <  i  <  j 
h  +  '<lj<Pj  (2-13) 

aj  <  li  <  Pi  for  t  >  j}. 
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A  A 

Note  that  Fj  =  0  if  J  <  m  and  lj  =  o \j  or  if  j  >  m  and  lj  =  Pj. 

What  remains  is  partitioned  into  disjoint  undetermined  intervals  Uj  given  by 


Uj  =  {/  £  U  :  li  <  U  <  pi  for  i  <  j 


h  <  h  <  h  - 1 

U  <U  <  Pi  for  j  <  i  <  m 

a,  <  li  <  /i  for  i  >  m},  if  j  <  m 


=  {!€(/: 


h  <  h  <  Pi  for  i  <  m 
<  li  <  li  for  m  <  i  <  j 


Ij  +  l<lj<  h 

a,-  <  li  <  h  for  i  >  j},  if  j  >  m. 


(2.14) 


Note  that,  as  before,  Uj  is  empty  if  lj  =  lj.  Otherwise,  Uj  will  be  maintained  for 
future  iterations.  (Again,  if  limiting  indices  were  not  found,  we  replace  /,  by  a,  if 
i  <m  and  by  $  if  i  >  m  in  the  description  of  Uj.) 

We  summarize  the  infeasible  decomposition  below  (again,  pre-screening  has  as¬ 
certained  that  the  interval  being  decomposed,  U  =  [a,/?],  contains  both  feasible  and 
infeasible  points). 


Procedure  INFEASIBLE (t/) 

Step  1  Start  with  the  state  ( p\a )  and,  while  S’(-)  is  not  satisfied,  determine  an 
infeasible  state  /  by  pushing  down  from  Pj  for  j  <  m  and  up  from  for 
j  >  m. 

Step  2  Determine  limiting  infeasible  indices  lj  as  in  Definition  2.9.  (If  feasible  sets 
are  not  to  be  removed,  set  i  =  Mfl.) 

Step  3  For  j  =  l,...,o  :  Form  Fj  as  in  (2.12)  or  (2.13).  If  Fj  ^  0,  set  P,F  = 
Pf  +  P{Fj),  where  P{Fj }  is  calculated  as  in  (2.2). 

25 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


Step  4  For  j  =  1,. . .  ,a  :  If  lj  ^  lj,  file  Uj  in  U,  where  U,  is  defined  by  (2.14). 

Finally,  we  give  the  entire  General  State  Space  Decomposition  algorithm. 

Procedure  GSD 
Step  1  Start  with  U  =  ft.  Set  U  =  0  and  PF  =  0. 

Step  2  Let  a  and  fi  denote,  respectively,  the  lower  and  upper  end  points  of  U.  If 
(a|/?)  is  not  feasible,  go  to  step  4.  If  (/?| a)  is  feasible,  set  P,F  =  PF  +  P{U} 
and  go  to  step  4. 

Step  3  Perform  either  FEASIBLE(IZ)  or  INFEASIBLE(tZ). 

Step  4  If  U  =  0,  stop.  Otherwise,  remove  a  set  U  from  U  and  go  to  step  2. 

The  statement  of  step  3  is  noticeably  vague.  How  should  we  decide  whether  to 
decompose  a  given  interval  feasibly  or  infeasibly?  We  consider  the  following  options: 

Pure  Option  1  All  intervals  are  decomposed  from  the  feasible  extreme. 

Pure  Option  2  All  intervals  are  decomposed  from  the  infeasible  extreme. 

Mixed  Option  1  Each  interval  is  decomposed  both  ways.  The  decomposition  for 
which  J2j=i  P(Uj)  is  smaller  is  used. 

Mixed  Option  2  Each  interval  is  decomposed  both  ways.  The  decomposition 
which  yields  fewer  new  undetermined  intervals  is  used. 

Mixed  Option  3  Each  interval  is  decomposed  both  ways.  The  decomposition  that 
produces  new  undetermined  intervals  containing  fewer  state  points  is  used. 

The  idea  behind  the  mixed  options  is  to  choose  the  “‘best”  decomposition  in 
each  iteration.  In  practice,  they  seem  to  perform  better  than  either  pure  strategy 
in  some  cases.  However,  it  also  seems  clear  that  trying  to  “optimize”  with  regard 


to  one  of  the  above  criteria  in  each  step  (filing  less  probability,  filing  fewer  points, 
filing  fewer  sets)  in  no  way  guarantees  the  most  efficient  possible  evaluation  of  PF . 
Some  or  all  of  the  five  options  will  be  compared  in  Chapter  3  for  minimum  spanning 
tree  problems  and  in  Chapter  4  for  minimum  cost  network  flow  problems. 

Another  practical  consideration  deals  with  the  handling  of  the  list  of  undeter¬ 
mined  sets.  Again,  a  number  of  possibilities  exist.  Alexopoulos  [1]  suggested  that 
LIFO  queue  descipline  outperforms  FIFO  in  keeping  the  list  small,  and  that  keeping 
a  list  sorted  by  set  probability  works  well  when  the  objective  is  the  computation  of 
tight  bounds.  We  have  had  very  good  results  keeping  a  sorted  list  initially  (with  the 
most  probable  undetermined  interval  on  top),  switching  to  LIFO  when  the  proba¬ 
bility  of  the  most  probable  filed  set  is  suitably  small.  All  numerical  results  in  +his 
thesis  use  this  form  of  list  management.  Sorting  the  list  by  probability  is  especia^y 
helpful  in  yielding  tight  bounds  on  PF  more  quickly,  as  discussed  in  Section  2.6. 

Finally,  we  have  stated  that  feasible  and  infeasible  intervals  produced  in  each 
iteration  are  discarded  at  the  end  of  that  iteration.  There  may  be  situations  in 
which  it  is  useful  to  save  one  or  both  of  these  set  types  in  files  for  future  use.  (Recall 
that  each  interval  is  stored  using  2a  integer  memory  locations  since  we  need  only  the 
lower  and  upper  end  points.)  For  example,  if  the  possible  values  Wj(-)  are  known 
for  each  X}  but  the  respective  probabilities  either  are  not  known  or  are  subject  to 
change,  we  may  retain  the  feasible  sets  (which  will  still  represent  a  partitioning  of 
the  feasible  region)  and  easily  re-calculate  PF  based  on  the  new  probabilities.  This 
is  very  helpful  in  performing  sensitivity  analyses.  At  other  times,  we  might  wish  to 
retain  the  infeasible  sets  for  further  analysis.  For  example,  in  Chapter  3  we  discuss 
passing  the  infeasible  sets  from  the  evaluation  of  the  probability  that  the  weight  of 
a  MST  does  not  exceed  d  as  input  to  the  routine  that  calculate  ci  ideality  indices 
so  as  to  compute  the  probability  that  an  arc  e  is  on  a  MST  whose  weight  exceeds  d. 
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2.5  Theoretical  Results  Concerning  GSD 

In  this  section,  we  digress  to  tie  up  some  theoretical  loose  ends  concerning  the 
General  State  Space  Decomposition  algorithm.  First,  we  demonstrate  the  connec¬ 
tion  between  our  algorithm  and  factoring.  Without  loss  of  generality,  we  consider 
decomposition  from  the  feasible  extreme.  In  each  iteration  we  decompose  an  un¬ 
determined  set  U  into  a  feasible  set  F,  infeasible  sets  sets  {Ij}j^A  and  new  unde¬ 
termined  sets  {f/j  -  Let  UF  denote  the  feasible  portion  of  U .  From  a  factoring 
standpoint,  then,  we  have  done  the  following. 

P{UF }  =  p{uFnF}  +  J2p{uFnij}  +  Y,p{uF  nt/j} 

i=i  J=1 

=  P{F}  +  o  4-  •  •  •  +  o  +  P{UF  n  [/,}  +  •••  +  P{UF  n  ua) 

=  P{F}  +  P{Uf}  +  ---  +  P{U[}. 

P{F }  is  added  to  the  accumulated  lower  bound  for  PF.  The  sets  Uj  are  filed  for 
future  iterations,  so  that  P{UF}  can  be  calculated  in  precisely  this  way. 

Secondly,  we  present  remarks  concerning  the  number  of  new  undetermined  inter¬ 
vals  produced  in  each  iteration  of  GSD.  In  particular,  we  shall  investigate  the  conse¬ 
quences  of  removing  various  combinations  of  feasible  and  infeasible  sets,  concluding 
that  the  methods  given  in  Section  2.4  cannot  be  improved  upon  in  this  regard.  We 
shall  look  chiefly  at  what  we  have  called  decomposition  from  the  feasible  extreme. 
By  symmetry,  the  results  also  apply  to  decomposition  from  the  infeasible  extreme. 
For  convenience,  and  without  loss  of  generality,  we  assume  in  this  section  that  all 
variables  are  lower  feasible,  and  that  we  seek  to  decompose  a  set  U  =  [a,  /?]  for 
which  or  is  a  feasible  point  and  fi  is  an  infeasible  point. 

Our  first  result  concerns  the  number  of  undetermined  intervals  needed  to  parti¬ 
tion  what  remains  of  U  after  the  removal  of  a  single  feasible  set  of  the  form  described 
in  Proposition  2.1.  The  symmetric  result  for  infeasible  sets  is  implied. 
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Proposition  2.2  Let  l  6  U  be  a  feasible  point  which  differs  from  the  point  0  in  i 

positions  (0  <  t  <  a),  and  let  F  be  the  feasible  set  defined  by  l  using  Proposition  2.1 . 

Then  exactly  i  intervals  are  required  to  partition  U\F. 

Proof:  We  have  already  seen  that  using  (2.10)  GSD  produces  one  undetermined 
interval  for  each  j  £  A  with  <  0r  In  the  present  setting,  then,  i  intervals  would 
be  produced.  Thus  we  need  only  show  that  it  is  not  possible  to  partition  U\F  using 
fewer  than  i  intervals.  We  shall  do  this  by  using  mathematical  induction  on  a,  the 
number  of  components. 

The  result  holds  for  a  =  2,  as  we  now  demonstrate.  If  i  =  0,  then  1  —  0  and 

U\F  =  0.  If  i  =  1,  then  l  differs  from  0  in  one  position,  and  hence  U\F  0, 

so  that  it  is  not  possible  to  cover  U\F  with  0  sets.  Finally,  if  i  =  2  then  clearly 
an  “L-shaped”  lattice  remains,  and  it  is  thus  impossible  to  cover  U\F  with  0  or  1 
interval. 

Now  suppose  the  result  holds  for  a  —  1  variables  (Xi, . . . ,  X„_i)  and  consider  the 
case  for  a  variables.  It  will  be  useful  to  view  0o,  the  state  space  of  (X,, . . .  ,  X0),  as 
consisting  of  n„  ordered  copies  of  the  state  space  of  (Xj, . . .  ,X0_j),  indexed 

by  (1, . . .  ,na},  and  to  view  undetermined  intervals  U  as  being  similarly  made  up  of 
ordered  lower-dimensional  undetermined  intervals. 

The  result  clearly  holds  when  i  =  0  or  i  =  1  by  precisely  the  same  logic  as  was 
used  above.  Now  suppose  that  2  <  t  <  a.  We  prove  that  it  is  impossible  to  partition 
U\F  using  fewer  than  i  intervals  by  contradiction.  In  light  of  our  previous  remark, 
U  will  be  viewed  as  consisting  of  ordered  (a  —  1  )-dimensional  undetermined  intervals 
designated  by  Uaa, . . .  ,Upa.  The  set  F  resides  in  sets  Uaa, . . .  ,Uja,  where  it  forms 
similar  sets  Faa,...  ,  Fjm  (where  by  similar  we  mean  that  each  Fj  is  oriented  with 
respect  to  the  other  points  in  Uj  in  the  same  way).  We  shall  consider  separately  the 
cases  /„  =  0a  and  la  <  0a. 

If  /„  =  0a,  then  clearly  *  <  o  — I  (/  agrees  with  0  in  at  least  one  position).  In  (//?„, 
the  point  J  differs  from  the  point  0  in  i  positions.  Any  partition  of  U\F  consisting  of 
fewer  than  i  intervals  clearly  implies  a  partition  of  Upa\Fpa  consisting  of  fewer  than 


i  intervals,  impossible  by  the  induction  hypothesis.  Thus  i  intervals  are  needed. 

If/.  <  /?.,  then  in  Uja  the  feasible  point  l  differs  from  the  point  (/?!,... , i,  /.) 
in  i  —  1  positions.  Now  suppose  that  there  is  a  partition  of  U\F  consisting  of  fewer 
than  i  intervals.  The  number  of  intervals  in  this  partition  must  be  i  —  1  since  a 
smaller  number  would  imply  an  impossible  partition  of  Uj a\Fja,  as  before.  fact, 
each  set  in  the  partition  must  intersect  Uja  yielding  a  non-empty  interval  that  does 
not  overlap  with  Fja.  Since  all  of  U\F  must  be  covered  by  this  partition,  some 
interval  must  contain  a  point  in  (Jja\Fja  and  the  point  (aj, . . . , a._i,  /?.),  which  lies 
“above”  F,  impossible  by  the  definition  of  an  interval.  So  i  intervals  are  needed  and 
the  result  is  proved.  □ 

An  important  special  case  of  Proposition  2.2  occurs  when  /  differs  from  the  point 
0  in  only  one  position,  in  which  case  U\F  is  an  interval.  Since  the  symmetric  result 
for  infeasible  sets  holds,  we  have  the  following  immediate  consequence. 

Corollary  2.1  Let  { Ij  }jc.a  be  a  collection  of  infeasible  sets,  each  of  which  is  defined 
by  an  infeasible  point  that  differs  from  the  point  a  in  only  one  position.  Turn 
U\(Uj^Ij)  is  an  interval. 

That  is,  as  long  as  we  restrict  ourselves  to  infeasible  sets  whose  defining  points 
differ  from  a  in  only  one  position,  we  may  remove  arbitrary  collections  of  infeasible 
sets  in  addition  to  the  feasible  set  F  without  impacting  the  number  of  undetermined 
intervals  produced.  Of  course,  a  symmetric  result  is  implied  for  removing  feasible 
sets  along  with  an  infeasible  set.  Thus  Corollary  2.1  clearly  applies  to  GSD  as  stated 
in  Section  2.4.  In  fact,  in  view  of  (2.10)  and  (2.14),  removing  F  from  the  reduced 
interval  //\(UJ€>/J)  will  likely  result  in  fewer  interval.  Thus  we  have  proved  the 
following  key  result. 

Theorem  2.1  In  GSD,  removing  infeasible  ( feasible )  sets  from  an  undetermined 
set  U  in  addition  to  the  primary  feasible  ( infeasible )  set  cannot  increase  the  number 
of  intervals  needed  to  partition  the  remainder  of  U . 
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Thus  GSD  is  able  to  remove  both  types  of  sets  simultaneously  without  splintering 
the  remaining  undetermined  points  into  many  sets.  One  might  wonder  whether  this 
fact  distinguishes  GSD  in  any  way,  or  if  indeed  any  routine  designed  to  remove 
feasible  and  infeasible  sets  would  perform  as  well  in  this  regard.  The  key  to  the 
answer  lies  in  Theorem  2.2  below.  (The  cases  in  which  iF  <  1  or  i1  <  1  are  already 
dealt  with  in  Preposition  2.2.) 

Theorem  2.2  Let  lF  be  a  feasible  point  which  differs  from  the  point  0  in  iF  >  1 
positions ,  and  let  F  be  the  corresponding  feasible  set.  Let  V  be  an  infeasible  point 
which  differs  from  the  point  a  in  i1  >  1  positions,  and  let  I  be  the  corresponding 
infeasible  set.  Then  the  number  of  sets  needed  to  partition  U\(F  U  I)  is  greater  than 
or  equal  to  the  number  of  sets  needed  to  partition  U\F. 

Proof:  (By  induction)  If  a  =  2,  then  we  clearly  have  iF  =  il  —  2.  Already 

Proposition  2.2  provides  that  U\F  needs  2  intervals  for  its  partition  (i.e.,  that  U\F 
is  an  “L-shaped”  lattice).  Since  i1  —  2  (i.e.,  V  lies  strictly  in  the  interior  of  U\F),  it 
is  clear  that  U\(F  U  7)  is  not  a  single  interval.  Thus  at  least  2  intervals  are  needed. 

Now  suppose  that  the  result  holds  for  a— l  variables  (X\, . . . ,  Xa_i ),  and  consider 
the  case  of  a  variables  (Xj,...  ,  X0).  As  before,  it  is  useCii  to  think  of  U  as  being 
composed  of  ordered,  indexed  (o  -  l)-dimensional  undetermined  sets  Ua,, . . . ,  Up,. 
We  consider  two  separate  cases. 

Case  1  ( lF,  >  lj.  for  some  j'  €  .4):  Without  loss  of  generality,  assume  that 
j*  =  a.  If  lF  =  0a,  then  also  Va  =  0a.  In  the  (a  —  1  )-dimensional  interval  Up,,  then, 
we  have  a  feasible  set  Fpa  =  F'l  Up,  whose  defining  point  lF  differs  from  0  in  iF 
positions,  and  an  infeasible  srt  Ipa  ~  In  Up,.  By  the  induction  hypothesis,  the 
partition  of  Up,\(Fp,  U  Ip,)  requires  at  least  iF  intervals.  Consequently,  at  least  iF 
are  needed  to  partition  the  set  U\(F  \J  I). 

Now,  suppose  that  lF  <  0a.  Ujf  contains  a  feasible  set  Fjf  —  F  0  U\f  whose 
defining  point  differs  from  {0\,. . .  ,0a~\  JF)  in  iF  -  1  positions  and  an  infeasible  set 
Ijf.  By  the  induction  hypothesis,  the  partition  of  Ujf\{FjfL)Iif)  requires  at  least  iF - 
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1  intervals.  Note  that,  by  the  definition  of  an  interval,  none  of  these  iF  —  1  intervals, 
even  when  extended  to  all  levels  of  U ,  can  contain  the  point  (orj , . . . ,  aa_i ,  If ),  which 
lies  “above”  the  set  F.  So  more  than  these  iF  —  1  intervals  are  needed  to  partition 
U\(F  U  I).  This  concludes  this  case. 

Case  2  (if  <  Jf  for  all  j  G  A):  In  this  case,  necessarily  If  <  f3j  for  all  j  G  A ,  so 
that  iF  =  a.  Similarly,  i1  =  a.  Suppose  that  there  exists  a  partition  of  U\(F  U  I) 
consisting  of  fewer  than  a  intervals.  This  partition  intersects  the  (a  —  l)-dimensional 
interval  Uaa  to  form  a  partition  of  Uaa\Faa  (where  Faa  =  F  D  Uaa)  consisting  of 
fewer  than  a  intervals.  Now  in  Uaa,  the  defining  point  of  Faa  differs  from  the  point 
(/?!,. ..  ,/30_i,a0)  in  a  —  1  positions,  and  hence  by  Proposition  2.2  a  —  1  intervals 
are  required  to  partition  Uaa\Faa.  Thus  the  partition  of  U\(F  U  / )  requires  a  —  1 
intervals,  each  of  which  intersects  Uaa  in  a  nonempty  interval  outside  of  the  set  F . 
This  means  that  one  of  these  intervals  must  contain  points  in  Uaa\Faa  and  the  point 
(aj , . . . ,  a0_ j ,  If  +  1 )  (a  point  “above”  F).  This  is  not  possible  given  the  definition 
of  an  interval.  Thus  at  least  iF  intervals  are  needed  to  partition  U\(F  U  /).  This 
concludes  this  case  and  the  result  is  now  proved.  □ 

Thus  we  see  that  removing  even  a  single  infeasible  set  not  of  the  form  used 
in  GSD  can  result  in  an  increase  in  the  number  of  sets  needed  to  partition  the 
remainder.  This  situation  would  only  be  compounded  by  the  removal  of  multiple 
arbitrary  infeasible  sets.  The  following  result  gives  an  upper  bound  on  the  number 
of  intervals  needed  to  partition  U\(F  U  I)  in  Theorem  2.2. 

Proposition  2.3  Let  F  and  I  be  as  in  Theorem  2.2.  Then  U\(F  U  I)  can  be 
partitioned  into  iF  +  i1  —  1  intervals. 

Proof:  Let  j'  G  A  be  such  that  If.  <  if.  (We  must  be  able  to  find  j *  since  F  and 
I  cannot  overlap.)  Define 

Ur  =  {/  G  U  :  If.  +  1  </>•</;.-  1 

<*i<li<0i  for  i  /  y*  } 


32 


(Note  that  Uj •  =  0  if  JJ.  =  lj,  —  1.).  What  remains,  then,  is  an  interval 

Uf  =  [<*,  {0\,  •  •  •  0j-+\,  ■  •  •  >/?a)] 

containing  the  feasible  set  F  (whose  defining  point  differs  from  UF's  upper  endpoint 
in  iF  —  1  positions),  and  an  interval 

Uj  =  [(ati, . . .  . . .  ,aa),^] 

containing  the  infeasible  set  /  (whose  defining  point  differs  from  Uhs  lower  endpoint 
in  i1  —  1  positions).  By  Proposition  2.2,  we  partition  Up  using  iF  —  1  sets  and  Uj 
using  i1  —  1  sets.  Thus,  counting  f/j.,  we  have  used  (iF  —  l)  +  (i/  —  1)  + 1  =  iF  +  il  —  1 
intervals  in  the  decomposition.  □ 

To  summarize,  our  goal  in  this  section  was  to  investigate  the  consequences  of 
removing  infeasible  sets  in  addition  to  a  primary  feasible  set  (or  feasible  sets  in 
addition  to  a  primary  infeasible  set).  In  GSD,  such  an  action  causes  the  number 
of  intervals  needed  for  the  partition  either  to  stay  the  same  or  to  decrease.  If  the 
secondary  sets  are  not  of  the  form  used  in  GSD  (i.e.,  if  they  are  infeasible  sets  whose 
defining  points  differ  from  a  in  more  than  one  position  or  if  they  are  feasible  sets 
whose  defining  points  differ  from  0  in  more  than  one  position),  the  number  needed 
either  stays  the  same  or  increases.  This  brings  us  to  the  following  result. 

Corollary  2.2  Let  V  be  the  class  of  all  decomposition  schemes  that  consist  of  re¬ 
moving  a  feasible  ( infeasible )  interval  and  one  or  more  infeasible  ( feasible )  intervals 
from  an  undetermined  interval  U  using  Proposition  2.1.  Then  GSD  is  minimal  in 
V  with  respect  to  the  number  of  new  undetermined  sets  produced  per  iteration. 

One  final  point  must  be  made.  Each  decomposition  scheme  in  V  involves  a 
variety  of  complex  interactions  that  determine  its  ultimate  effectiveness  on  a  given 
problem.  Experience  indicates  that  no  single  objective  (e.g.,  keeping  the  list  of 
sets  short,  removing  maximum  probability,  expending  little  computational  effort) 
can  guarantee  across-the-board  success.  Trying  to  file  few  sets  does  seem  to  play  a 
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big  role,  and  intuitively  contributes  in  avoiding  the  splintering  of  a  set  into  many 
small  subsets.  Nevertheless,  because  of  the  complex  interplay  mentioned  above,  we 
shall  attempt  in  Chapters  3  and  4  to  remove  a  single  arbitrary  infeasible  set  and  an 
arbitrary  feasible  set  in  each  iteration,  ignoring  the  warning  of  Corollary  2.2.  This 
“hybrid”  approach  will  be  compared  with  pure  and  mixed  GSD  applications  for  the 
stochastic  minimum  spanning  tree  and  minimum  cost  flow  problems. 

In  the  final  sections  of  this  chapter,  we  conclude  our  description  of  GSD  by 
detailing  its  use  for  producing  bounds  and  providing  input  for  variance-reducing 
Monte  Carlo  estimation  schemes. 

2.6  Using  GSD  to  Produce  Bounds 

In  Section  2.4  we  presented  an  algorithm  for  computing  PF  exactly  in  a  number 
of  iterations  that  is  typically  small  (relative  to  |fl|)  for  moderate  size  problems.  For 
large  problems,  this  may  still  represent  an  unacceptable  computational  workload. 
In  such  cases  the  real  strength  of  a  routine  like  GSD  lies  in  its  ability  to  generate 
information  after  each  iteration  which  can  be  used  for  computing  bounds,  design¬ 
ing  efficient  Monte  Carlo  sampling  plans,  and  conducting  sensitivity  analysis  with 
respect  to  changing  input  distributions  [2], [4].  We  have  pointed  out  that  GSD  has 
available  at  all  times  a  lower  bound  for  PF .  In  this  section,  we  briefly  discuss  how 
we  can  stop  short  of  complete  decomposition  of  f!  and  retrieve  both  lower  and  upper 
bounds  for  PF . 

Suppose  we  carry  out  a  number  of  iterations  of  GSD,  resulting  in  accumulated 
feasible  probability  PF  and  a  collection  of  m  remaining  disjoint  undetermined  in¬ 
tervals  which  we  may  renumber  as  {C/<  Given  the  way  it  was  accumulated,  it  is 
clear  that  PF  is  a  lower  bound  for  PF .  Since  all  feasible  points  not  yet  accounted 
for  are  contained  in  U”1  ,(/<,  we  have 

m 

pf  <  pf  <  pf  +  Y.p{Vi}  =  pu-  (2.15) 

i=i 
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That  is,  stopping  the  GSD  procedure  at  any  point  yields  lower  and  upper  bounds  on 
the  probability  of  interest.  Furthermore,  since  nearly  every  iteration  removes  some 
undetermined  probability  (by  identifying  feasible  and  infeasible  sets)  and  adds  to 
the  accumulated  feasible  probability,  the  sequence  of  lower  bounds  is  nondecreasing 
and  the  sequence  of  upper  bounds  is  nonincreasing.  The  decomposition  may  then 
be  terminated  based  on  one  of  the  following  stopping  conditions: 

•  Stop  when  a  desired  number  of  sets  has  been  decomposed. 

•  Stop  when  bounds  are  sufficiently  close. 

We  may  alter  the  procedure  GSD  above  by  substituting  step  4  for  step  4  and  adding 
steps  5  and  6: 

Step  4  If  U  =  0,  stop.  If  a  desired  stopping  condition  is  met,  go  to  step  5.  Other¬ 
wise,  remove  a  set  U  from  U  and  go  to  step  2. 

Step  5  Let  PF  =  PF . 

Step  6  For  i  =  1  Evaluate  extremes  of  f/<.  If  Ui  is  infeasible,  discard  U ,. 

If  Ui  is  feasible,  set  Pf  =  Pf  +  P{U,}  and  Pf  =  Pf  +  P{Ut}.  Otherwise, 
since  Ui  contains  some  feasible  points  but  is  not  completely  feasible,  set  Pf  = 

pf  +  pm- 

(In  step  6,  we  may  skip  the  evaluation  of  the  extremes  and  simply  set  Pf  = 
Pf  +  P{U<}  evaluating  the  extremes  is  deemed  computationally  expensive.)  After 
the  revised  GSD  has  been  performed,  we  are  left  with  lower  bound  Pf  and  upper 
bound  Pf,  and  with  a  set  of  m  <  m  undetermined  intervals.  These  bounds  and 
remaining  undetermined  sets  may  be  used  as  input  to  Monte  Carlo  sampling  plans 
for  estimating  PF,  as  described  in  [1],  [2j  and  presented  in  the  following  section. 

The  revised  GSD  is  important  since  in  many  cases,  deriving  a  set  of  tight  bounds 
is  just  as  useful  as  evaluating  PF  exactly  (especially  for  large  problems,  in  which  it 
is  the  only  practical  recourse).  In  fact,  it  is  the  desire  for  quick,  tight  bounds  that 
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motivated  some  of  the  features  of  GSD.  For  example,  filing  and  retreiving  sets  by 
probability  is  an  attempt  to  ensure  that  no  high-probability  sets  remain  undeter¬ 
mined  (an  important  goal  since  bounds  can  differ  by  as  much  as  P{Ui})-  Some 
decomposition  approaches  do  not  identify  limiting  indices  lj  so  that  only  a  feasible 
set  and  undetermined  sets  are  produced  in  each  iteration.  Thus  the  upper  bound  is 
reduced  only  by  intervals  which  are  found  to  be  completely  infeasible.  When  lim¬ 
iting  indices  are  used,  both  bounds  may  improve  in  nearly  every  iteration.  Finally, 
no  existing  decomposition  approach  decomposes  both  feasibly  and  infeasibly.  We 
anticipate  that  for  some  problems  this  added  flexibility  will  allow  as  much  tightening 
of  bounds  as  possible  from  both  directions,  with  benefits  exceeding  the  cost  of  the 
additional  computational  effort. 

2.7  GSD  Bounds  as  Input  to  Monte  Carlo  Esti¬ 
mation  Routines 

We  may  take  the  partial  decomposition  results  of  Section  2.6  and  use  them 
as  input  to  a  highly  efficient  Monte  Carlo  sampling  routine.  The  sampling  plan 
described  in  this  section  combines  importance  and  stratified  sampling,  requires  no 
more  work  per  iteration  than  a  crude  Monte  Carlo  sampling  of  fl,  and  has  much 
smaller  variance  than  such  a  plan.  In  fact,  the  ratio  of  crude  Monte  Carlo  variance 
to  the  variance  of  this  sampling  plan  based  on  GSD  output  is  at  least 


1/  f(i-r,n]\ 


(2.16) 


which  can  be  computed  before  sampling  begins. 

We  begin  with  bounds  Pf  and  Pf,  and  m  remaining  undetermined  intervals 
{f/,}^!.  Let  7T,  =  P{Ui}  for  all  i  and  let  it  =  £”L,  Wi.  For  /  €  0,  define  <f>(l)  =  I  if  / 
is  a  feasible  point,  <f>(l)  =  0  if  /  is  infeasible.  Suppose  we  wish  to  draw  N  independent 
samples  from  the  undetermined  sets.  For  i  =  l,...,rh  we  shall  draw  N,  =  Nni/n 
samples  from  set  U,  as  follows. 
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For  n  =  l,..., N^. 

•  For  each  j  =  1, . .  .  ,a:  sample  a  level  l}  according  to  the  distribution 


» 


ft 


•  Form  the  resulting  point  =  (li, . . . ,  la)  and  evaluate 
Then 

m  Ni 

PF  =  Pf  +  £>•  £  ^(X<“>)/JV,  (2.18) 

i=l  |n=l 

is  an  unbiased  estimator  of  PF  with  variance  as  given  in  [1]. 

This  efficient  sampling  routine  is  easy  to  implement,  as  the  U,'s  may  be  taken 
one  at  a  time,  used  for  sampling  and  then  discarded.  It  adds  even  greater  flexibility 
to  the  GSD  process,  especially  for  large  problems. 

2.8  Conclusions 

In  this  chapter,  we  have  presented  the  General  State  Space  Decomposition  algo¬ 
rithm,  a  powerful  and  flexible  tool  for  evaluating  measures  related  to  the  operation 
of  stochastic  systems.  Starting  with  the  significant  groundwork  laid  by  previous 
studies,  GSD  has  strengthened  the  original  decomposition  of  Doulliez  and  Jamoulle 
in  the  following  ways: 

1)  The  general  framework  and  vocabulary  developed  in  Section  2.2  give  a  common 
point  of  departure  for  discussing  and  comparing  the  application  of  these  types 
of  decompositions.  The  GSD  routine  itself  represents  an  advancement  as  an 
effective  application  of  factoring  that  is  general  enough  to  be  easily  tailored 
to  address  a  variety  of  stochastic  problems. 
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The  content  of  Section  2.5  is  important  in  establishing  GSD  as  best  (in  min¬ 
imizing  the  number  of  undetermined  sets  produced  in  each  iteration)  among 
similar  alternative  approaches. 

Previous  decomposition  approaches  employed  either  Pure  Strategy  1  or  Pure 
Strategy  2  (see  Section  2.4).  Such  strategies  will  likely  be  good  for  some  choices 
of  constraints,  bad  for  others  (see  Section  3.4).  GSD  as  applied  in  Chapter 
3  represents  the  first  routine  of  this  type  which  allows  mixed  strategies,  thus 
attempting  to  gain  the  benefits  of  both  feasible  and  infeasible  decomposition 
when  possible. 

The  use  of  GSD  both  in  Chapter  3  and  in  Chapter  4  continues  to  demonstrate 
the  importance  of  efficient  list  processing  for  U ,  being  one  of  the  first  such 
applications  to  keep  the  undetermined  sets  in  a  sorted  heap  structure  (see  also 


CHAPTER  3 


The  Stochastic  Minimum 
Spanning  Tree  Problem 


3.1  Introduction 

Consider  an  undirected,  connected  graph  G  =  (Af,A)  with  node  set  Af  = 
{l,...,n}  and  arc  set  A  =  {l,...,a}.  A  spanning  tret  is  a  maximally  acyclic 
connected  subgraph  of  G,  so  named  because  it  reaches  every  node  in  Af.  Suppose 
that  each  arc  in  A  has  an  associated  nonnegative  weight.  An  important  problem 
in  network  optimization  is  the  identification  of  a  spanning  tree  with  the  smallest 
possible  total  weight.  This  so-called  minimum  spanning  tree  (MST)  problem  is  im¬ 
portant  in  its  own  right  as  it  attempts  to  connect  all  nodes  in,  e.g.,  a  communications 
network  as  economically  as  possible.  In  addition,  there  are  other  unrelated  network 
optimization  problems  that  can  be  cast  as  MST  problems.  A  number  of  efficient 
algorithms  exist  for  solving  the  minimum  spanning  tree  problem  (see  (24], [28], [29]). 

In  this  chapter,  we  are  concerned  with  describing  the  probabilistic  behavior  of 
MSTs  in  networks  whose  arcs  have  random  finite  weights.  Specifically,  the  weight 
of  arc  j  £  A  is  a  discrete  random  variable  Xj  that  takes  on  weight  values  u>j(l)  < 
...  <  Wj(nj)  with  respective  probabilities  Pj(l), . . .  ,Pj(rij).  We  assume  that  the 
components  of  X  =  (Xi, . . .  ,Xa)  are  mutually  independent.  The  state  space  D 
of  X  contains  |ftj  =  f]“=i  nj  elements  of  the  form  x  =  (u>i (/i wa(la)),  where 
lj  €  {1,. . .  ,rij}  for  each  j.  As  in  Section  2.2,  we  employ  the  convention  of  writing 
this  state  point  as  /„). 
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We  shall  concern  ourselves  in  this  chapter  with  the  following  problems  related 
to  MSTs  in  stochastic  networks. 

Problem  STl  Define  the  random  variable  W(X)  as  the  weight  of  a  MST  and  let 
d  >  0.  We  want  to  compute  the  probability  P{VF(X)  <  d }  that  we  can 
construct  a  spanning  tree  whose  weight  does  not  exceed  d. 

For  this  problem  a  realization  /  of  X  is  considered  feasible  if  W(l)  <  d,  and  clearly 
Xi,...,Xa  are  all  lo.ver  feasible.  This  problem  will  be  discussed  in  Section  3.2, 
including  a  modification  to  GSD  that  will  allow  the  simultaneous  determination  of 
the  entire  distribution  of  l'F(X). 

Problem  ST2  Let  e  G  A.  We  want  to  compute  the  probability  that  e  belongs  to 
a  minimum  spanning  tree. 

Associated  with  every  realization  /  of  X  is  a  collection  T(l)  =  {  7, , . . . ,  7'm  } 
of  minimum  spanning  trees.  The  state  point  l  is  feasible  if  e  G  Tx  for  some  i  G 
For  this  problem,  the  random  variable  Xe  is  lower  feasible  while  the 
weight  of  all  other  arcs  is  upper  feasible.  This  problem,  along  with  two  modifications, 
will  be  discussed  in  Section  3.3. 

Before  solving  these  problems  by  using  GSD,  it  will  be  instructive  to  make  some 
observations  about  the  dynamics  of  MSTs  when  an  arc  changes  weight.  These 
observations  will  be  useful  since  each  GSD  procedure  will  consist  of  a  sequence  of 
steps  in  which  the  level  of  a  single  arc  is  changed.  By  knowing  in  advance  the 
consequences  of  these  steps,  we  can  maintain  or  achieve  properties  such  as  tree 
membership  or  tree  weight. 

We  begin  at  an  arbitrary  point  l  =  (/lt...,/0).  Let  T  be  a  MST  at  this  level. 
We  shall  independently  decrease  and  increase  the  weight  of  arcs  both  on  and  off  the 
tree  T  addressing  the  following  questions  for  each  action: 

1)  Does  this  action  cause  T  to  cease  being  a  MST,  thus  forcing  a  change  in  tree 
membership? 
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2)  Does  this  action  change  the  weight  of  a  MST? 

We  summarize  these  observations  below  as  Properties  1-4.  For  convenience,  we 
shall  use  a  shorthand  notation  for  some  set  operations.  For  example,  the  operation 
T  U  will  be  written  as  T  -f  i  —  j. 

Property  1  (Decreasing  the  weight  of  an  arc  j  €  T)  Since  T  is  a  MST  with 
the  current  arc  weights,  it  will  surely  remain  so  if  the  weight  of  one  of  its 
arcs  decreases.  Of  course,  the  weight  of  T  will  decrease  by  the  amount  of  the 
reduction  in  the  weight  of  arc  j. 

Property  2  (Decreasing  the  weight  of  an  arc  j  ^  T)  If  j  £  T,  then  decreas¬ 
ing  its  weight  sufficiently  can  cause  it  to  displace  from  T  an  arc  with  greater 
weight.  Thus  we  must  first  investigate  C(j),  the  unique  cycle  that  arc  j  forms 
with  T.  Only  arc  i,  the  most  heavily  weighted  arc  in  C(j)  —  j,  can  be  so 
displaced  in  order  to  obtain  a  minimum  spanning  tree.  As  long  as  the  weight 
of  arc  j  remains  at  least  as  large  as  that  of  arc  i,  the  reduction  of  j  impacts 
neither  T  nor  its  weight.  Otherwise,  T  —  i  +  j  is  the  new  MST,  with  weight 
less  than  that  of  T. 

Property  3  (Increasing  the  weight  of  an  arc  j  T)  Since  j  £  T  with  current 
weights,  the  weight  of  arc  j  can  be  increased  arbitrarily  without  impacting  T 
or  its  weight. 

Property  4  (Increasing  the  weight  of  an  arc  j  £  T)  Increasing  the  weight  of 
j  can  cause  T  to  cease  being  a  MST.  Among  those  arcs  not  on  T  whose  cycles 
with  T  contain  j ,  let  i  be  the  arc  with  the  smallest  weight.  As  long  as  the 
weight  of  arc  j  does  not  exceed  that  of  arc  i,  T  is  still  a  MST  (though  with 
increased  weight).  Otherwise,  T  —  j  +i  is  a  MST  and  further  increases  in  the 
weight  of  j  are  governed  by  Property  3. 

We  shall  often  refer  to  the  process  of  exchanging  arcs  between  a  MST  and  the 
remainder  of  A  as  pivoting. 
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We  now  proceed  to  the  solution  of  Problems  ST1  and  ST2  in  Sections  3.2  and 
3.3,  respectively.  Section  3.4  contains  numerical  examples.  In  Section  3.5  we  suggest 
extensions  of  the  MST  results  to  other  matroid  settings.  Concluding  remarks  are 
made  in  Section  3.6. 

3.2  Computing  the  Probability  Distribution  of 
the  Weight  of  a  MST 

We  consider  here  the  stochastic  network  system  that  operates  feasibly  if  the 
constraint  VU(X)  <  d  is  satisfied,  and  seek  to  compute  the  probability  P{!U(X)  < 
d }.  Unfortunately  this  problem,  like  many  that  fit  the  general  problem  description 
in  Section  2.3,  is  provably  hard.  We  show  this  by  reduction  from  the  following 
#P-complete  problem  (see  [6]). 

All-terminal  Undirected  Rational  Network  Reliability  Problem  (ATRP) 

Inputs;  Undirected  network  with  arcs  that  function  or  fail  independently  of  each 
other.  For  arc  i ,  the  probability  of  failure  is  bi/ci,  where  b,  is  a  non-negative 
integer  and  c,  is  a  positive  integer  such  that  6,  <  c*. 

Outputs:  6  and  c,  where  c  is  a  positive  integer,  6  is  a  non-negative  integer  such 
that  b  <  c,  and  P{the  system  functions}  =  P{the  network  is  connected} 

=  P{there  is  a  spanning  tree  consisting  of  functioning  arcs  }  =  b/c. 

Proposition  3.1  The  evaluation  of  P{W(X)  <  d}  is  #P-hard. 

Proof:  Consider  an  instance  of  ATRP.  To  each  arc  i,  assign  the  random  weight 

I  0  with  probability  1  — 

1  (  1  with  probability  ^ 
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where  nonzero  weight  corresponds  to  failure.  Then  evaluating  P-fW^X)  <  0}  is 
equivalent  to  finding  the  probability  that  there  is  a  spanning  tree  made  up  of  func¬ 
tioning  arcs.  That  is,  P{W(X.)  <  0}  is  equal  to  the  all-terminal  network  reliability. 
Thus,  the  ability  to  evaluate  P{W(X.)  <  d }  implies  the  ability  to  solve  ATRP.  Since 
ATRP  is  #P-complete,  the  evaluation  of  P{W(X.)  <  d}  is  #P-hard.  □ 

As  discussed  earlier,  all  arc  weight  random  variables  are  lower  feasible.  The 
extreme  feasible  candidate  of  an  undetermined  set  U  =  [a,  /?]  is  simply  the  point  a 
and  the  extreme  infeasible  candidate  is  /?. 

We  are  now  in  a  position  to  discuss  the  application  of  GSD  to  the  present  prob¬ 
lem.  Given  the  groundwork  laid  in  Chapter  2,  this  will  amount  to  little  more  than 
specifying  the  mechanics  of  finding  l  and  the  1,' s  for  feasible  and  infeasible  decom¬ 
position.  As  remarked  in  Chapter  2,  there  are  several  ways  to  describe  feasible 
and  infeasible  decomposition  for  this  problem.  In  Section  3.2.1  we  present  the  end- 
product  of  a  long  evolution  of  decomposition  schemes.  It  is  fairly  involved,  but  has 
shown  itself  to  be  a  flexible  routine.  Section  3.2.2  contains  an  alternative  algorithm 
which  is  conceptually  much  simpler,  and  which  out-performs  the  more  sophisticated 
algorithm  of  Section  3.2.1  in  some  instances.  Section  3.2.3  contains  a  “hybrid”  rou¬ 
tine  based  on  the  decomposition  discussed  in  the  proof  of  Proposition  2.3.  Finally, 
Section  3.2.4  presents  a  modification  to  the  work  in  Section  3.2.1  that  allows  us  to 
compute  the  entire  distribution  of  W(X)  simultaneously. 

3.2.1  Spanning  'free  Decomposition  Algorithm  I 

The  algorithm  presented  in  this  section  is  designed  to  carry  out  both  feasible 
and  infeasible  decomposition,  as  outlined  in  Chapter  2.  It  can  be  executed  in  several 
ways,  namely  using  Pure  Strategies  1  and  2  and  Mixed  Strategies  1  through  3,  as 
described  in  Section  2.4.  All  five  options  will  be  used  to  solve  the  sample  problems 
of  Section  3.4. 

We  now  give  the  mechanics  of  Spanning  Tree  Decomposition  Algorithm  I  (STDA- 
I)  feasible  decomposition  and  of  STD  A- 1  infeasible  decomposition.  In  both  cases  we 


describe  a  single  iteration,  in  which  we  partition  an  undetermined  interval  U  — 
[a,  /?],  which  contains  feasible  and  infeasible  states. 


STDA-I  Feasible  Decomposition 

To  determine  the  feasible  point  /,  we  begin  at  the  extreme  feasible  point  a  and 
increase  arc  weights  while  remaining  feasible.  For  all  j  G  A,  we  set  the  arc  weight 
at  Wj(aj)  and  find  a  MST,  which  we  denote  T(a).  It  is  clear  by  Property  3  that  any 
j  not  in  T(a)  can  be  raised  to  weight  Wj{R3)  without  impacting  T(a).  Thus  with 
no  additional  work  we  may  set 


ft  if  j*T(a) 
Oj  if  j  G  T(a ) 


yielding  a  feasible  point  for  which  a  —  n  -f  1  arcs  .  >t  as  deeply  as  possible.  At 
this  point  we  could  try  to  increase  the  weight  of  some  of  the  arcs  on  the  tree  T(a) 
using  Property  4.  Experimentation  has  indicated  that  the  considerable  extra  work 
that  would  be  involved  does  not  pay  off  in  this  case.  Therefore,  we  shall  use  /  as 
given  above. 

All  that  remains,  then,  is  to  determine  the  limiting  feasible  indices  lj.  We  already 
have  determined  a  —  n  +  1  of  them:  since  h  >  for  all  j ,  our  assignment  of  / 
above  indicates  that  lj  =  for  all  j  ^  T(a).  Increasing  the  weight  of  the  tree  arcs 
requires,  according  to  Property  4,  some  pr  processing  to  determine  what  restrictions 
there  may  be  on  such  an  action.  Fortunately,  though,  since  each  /,  is  determined 
independently  and  starting  at  the  same  poi.  *  a,  all  the  tree  arcs  can  be  preprocessed 
in  one  pass  through  the  arcs  off  T( a)  by  the  following  procedure. 


Procedure  SHORT 


Step  1  For  j  G  T(a):  SHORT(j)  =  oc. 

Step  2  For  j  0  ^(o):  Find  C(j),  the  unique  cycle  j  forms  with  T(n).  For  all 
i  G  C(j)  — j ,  put  SHORT(i)  =  min{SHORT(i),w}{c»})}. 


Step  3  End. 


For  all  j  €  T(a)  this  procedure  produces  SHORT(j),  the  weight  of  the  off- 
tree  arc  that  would  displace  arc  j  if  the  weight  of  this  arc  were  to  increase  beyond 
S HO RT(j).  For  j  $  T(a)  we  set  lj  =  0}.  We  then  proreed  to  find  independently 
for  each  j  €  T(a)  as  follows.  We  begin  with  the  MST  T(a )  of  weight  W(a)  and 
all  arcs  at  levels  a  (this  information  having  been  saved  from  the  determination  of 
/).  If  lj  =  flj,  we  set  lj  =  flj  and  go  to  the  next  arc.  Otherwise,  we  set  lj  — 
max{/  :aj<l<0j ,  Wj(l)  <  SHORT(j)}. 

The  determination  of  /  involves  relatively  little  work.  Procedure  SHORT  takes 
time  O(an),  and  the  subsequent  determination  of  the  s  (which  is  done  for  the  tree 
arcs  only)  takes  time  0(n  *  maxjg^n^). 

It  might  seem  inconsistent  to  the  reader  that  in  finding  l  we  declined  to  raise 
any  tree  arcs  so  as  to  avoid  costly  preprocessing,  yet  we  did  the  preprocessing  for 
the  lj' s.  The  reason  for  this  apparent  inconsistency  is  that  these  situations  differ  in 
two  important  ways.  First,  with  the  lj' s  one  preprocessing  pass  through  the  off-tree 
arcs  was  sufficient  to  set  up  all  lj  determinations  because  the  weight  increases  used 
to  find  the  lj' s  were  independent  and  started  at  the  same  point  a.  With  /,  since 
the  point  is  determined  considering  all  arcs  jointly,  none  of  the  weight  increases  are 
independent.  And  since  in  general  each  arc  is  raised  from  a  different  state  point 
(wherever  the  previous  arc  left  off),  we  would  have  had  to  execute  SHORT  as  many 
as  n  —  2  times  (the  preprocessing  for  the  first  tree  arc  raised  was  done  in  finding  /). 
Secondly,  as  mentioned  in  Section  2.4,  lj  must  be  pushed  as  far  as  possible,  so  that 
there  was  no  choice  but  to  preprocess. 

After  limiting  feasible  indices  have  been  determined  for  all  arcs,  we  are  ready  to 
decompose  as  in  Section  2.4  and  proceed  to  the  next  undetermined  interval. 

STDA-I  Infeasible  Decomposition 

Infeasible  decomposition  begins  at  the  extreme  point  0  (which  by  pre-screening 
wc  know  to  be  an  infeasible  point)  and  pushes  in,  identifying  an  infeasible  point  / 
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and  limiting  infeasible  indices  lj.  As  before,  the  determination  of  for  each  j  and 
the  determination  of  /,  each  beginning  at  0,  can  be  designed  to  be  carried  out  in 
any  order.  In  feasible  decomposition,  we  found  /  first.  Here,  the  limiting  infeasible 
indices  lj  will  be  found  first.  Afterward,  we  use  the  fact  that  h  <  l  I  for  j  £  A  during 
the  derivation  of  /. 

We  begin  by  finding  T(0),  a  MST  at  the  extreme  infeasible  point,  with  weight 
W(0).  For  j  £  T(0)  we  can  use  Property  1  to  find  lj  with  relative  ease,  since  we  can 
look  at  arc  j  in  isolation,  ignoring  all  othei  arcs.  We  initially  set  l:  =  / 3j  and  lower 
one  level  at  a  time  until  either  the  MST  weight  would  drop  below  d  with  further 
reduction  or  lj  hits  cij  while  the  tree  weight  is  still  infeasible. 

On  the  other  hand,  finding  lj  for  off-tree  arcs  j  is  complicated  by  the  fact  that 
reducing  the  weight  of  arc  j  could  force  it  into  the  MST,  as  discussed  above  in 
Property  2.  So  we  identify  C(j),  the  unique  cycle  that  arc  j  forms  with  T(0), 
and  identify  i,  the  most  heavily  weighted  arc  (at  level  0)  in  C(j)  —  j.  We  then 
decrease  the  weight  of  arc  j  one  level  at  a  time  until  either  we  hit  and  still  have 
tWj(aj)  >  Wi(0i),  in  which  case  l}  =  aj  (here  we  also  preset  lj  =  otj ) ,  or  at  some 
point  the  weight  of  arc  j  becomes  strictly  less  than  that  of  arc  i  and  T{0)  —  i+j  is 
the  new  MST  (here  we  preset  lj  at  the  index  prior  to  that  which  forced  j  to  enter 
the  tree).  If  this  tree  weight  is  still  feasible,  we  continue  reducing  the  weight  of  arc 
j  as  for  tree  arcs.  If  not,  then  lj  is  set  at  the  level  just  prior  to  the  level  at  which 
arc  j  entered  the  tree. 

It  remains  to  derive  the  infeasible  point  /  by  reducing  weights  from  the  level  0. 
Rather  than  starting  from  1  =  0,  we  begin  with  the  point  given  by 

j  _  f  0,  'f  J  €  T(0) 

3  (  as  set  above  if  j  ^  T(0 ) 

It  is  clear,  by  Properties  1  through  4,  that  T(0)  is  a  MST  for  /  defined  in  this  way, 
and  that  W(0)  is  its  weight.  We  may  proceed  from  this  point  to  derive  a  deeper 
infeasible  point.  As  mentioned  earlier,  l}  >  lj  for  all  j  £  A.  Thus  the  state  /  serves 
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as  a  floor  for  the  further  reduction  of  l.  In  fact,  it  is  likely  that  for  some  arcs  j  we 
already  have  lj  =  lj  with  l  as  defined  by  (3.1).  Those  arcs  will  not  be  considered 
for  further  reduction.  The  remaining  arcs  will  be  reduced  one  at  a  time  as  much  as 
possible,  governed  by  the  applicable  Property,  as  long  as  the  tree  weight  does  not 
reach  d. 

This  further  reduction  could  be  accomplished  in  many  different  ways.  The  idea 
is  to  reduce  as  many  arcs  as  possible,  and  to  reduce  them  as  much  as  possible.  We 
want  to  avoid  using  arcs  j  (E  T(0)  for  which  Wj(lj)  —  Wj(lj  —  1)  is  large  since  such 
arcs  might  bring  the  MST  weight  down  to  d  all  at  once,  preventing  the  reduction  of 
other  arcs.  In  addition,  we  would  like  to  give  preference  to  arcs  with  a  Lt  of  room 
I]  —  lj  for  reduction.  In  order  to  try  to  accomodate  both  goals,  we  use  an  ordering 
of  arcs  that  is  inspired  by  a  knapsack  heuristic.  (If  n_,  =  2  for  all  j,  we  solve  a  0  —  1 
knapsack  problem.  The  space  taken  by  each  arc  in  that  case  is  Wj(/3j)  —  Wj(otj), 
all  arcs  have  equal  value,  and  the  available  space  in  the  knapsack  is  W(0)  —  d.) 
Specifically,  after  eliminating  arcs  with  lj  =  lj  from  consideration,  we  first  consider 
tree  arcs,  followed  by  arcs  not  on  T(/3),  in  decreasing  order  of  the  ratio  (lj  —  lj)  / 
[Wj(Jj)  —  Wj(lj  —  1)].  The  weight  of  each  arc  selected  for  reduction  according  to  this 
rule  is  decreased  one  level  at  a  time  until  further  reduction  either  is  not  possible 
(i.e.,  lj  =  lj)  or  would  result  in  an  infeasible  point. 

Once  /  has  been  determined,  we  are  ready  to  decompose  as  in  Section  2.4  and 
proceed  to  the  next  undetermined  interval. 

3.2.2  Spanning  Tree  Decomposition  Algorithm  II 

Section  3.2.1  gave  one  possible  application  of  GSD  to  the  determination  of  an 
ordinate  of  the  probability  distribution  of  minimum  spanning  tree  weights.  As  men¬ 
tioned  before,  there  is  nearly  limitless  flexibility  in  applying  GSD  to  a  given  problem 
instance.  STDA-I,  given  in  Section  3.2.1,  contains  what  seem  to  be  the  most  effec¬ 
tive  features  of  a  number  of  earlier  schemes.  This  section  describes  Spanning  Tree 
Decomposition  Algorithm  II  (STDA-II),  an  algorithm  that  is  conceptually  much 


47 


simpler.  We  include  it  because  it  is  fast  and  because  it  demonstrates  the  flexibil¬ 
ity  of  GSD  by  its  departure  from  the  routine  described  in  Chapter  2.  Differences 
include  a  modified  prescreening  phase  and  a  simplified  decomposition  which  results 
from  not  finding  limiting  indices  lj.  Also,  except  for  undetermined  intervals  ruled 
infeasible  in  prescreening,  no  sets  will  be  decomposed  infeasibly. 

Before  describing  STDA-II  we  can  simplify  the  descriptions  given  earlier  that 
define  the  partitioning  of  an  undetermined  interval.  Since  no  variables  are  upper 
feasible  and  no  limiting  indices  are  found,  we  have  simply 

F  =  {/  €  U  :  atj  <  I3  <  lj  for  j  G  A)  (3.2) 

and 

Uj  =  {/  G  U  :  a,  <  l{  <  li  for  i  <  j  (3-3) 

/,  +  !</,<  ft 

a,  <  /,  <  A  for  i  >  j} . 

STDA-II  can  be  stated  as  follows,  with  explanation  following. 

Procedure  STDA-II 

Step  1  Start  with  U  =  fl.  Set  U  =  0,  P[  =  0. 

Step  2  Let  a  and  0  denote  the  lower  and  upper  endpoints  of  U ,  respectively.  Find 
T(a),  a  MST  at  a,  with  weight  W(a).  If  W(a)  >  d,  then  U  is  infeasible;  go 
to  step  7. 

Step  3  Let  Wp(ct)  =  J2jeT{a)wj(ft)-  If  IF p(cx )  <  d,  then  W(0)  <  d  and  U  is 
feasible.  Set  Pf  =  P[  +  P{U)  and  go  to  step  7. 

Step  4  Set  lj  =  ft  for  j  $  T(a),  l3  =  a,  for  j  <E  T(a),  STWT  =  W(a).  Sort  the 
arcs  j  G  T(a)  for  which  atj  <  0}  in  increasing  order  of  Wj(0j)  — 
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Step  5  Taking  the  tree  arcs  in  this  order,  do  for  j  chosen:  Set  lj  =  0j,  STWT  = 
STWT  -f  Wj(Pj)  —  Wj(atj).  As  long  as  STWT  <  d,  continue  choosing  arcs. 
When  arc  j  causes  STWT  >  d,  put 

lj  =  max{7  :  atj  <7  <  0j,  STWT  -  Wj(0j)  +  Wj( 7)  <  d }  (3.4) 

and  go  to  step  6. 

Step  6  Form  F  as  in  (3.2)  and  set  Pf  —  Pf  +  P{F}.  For  j  £  A  :  If  lj  ^  0j,  file 
Uj  in  U,  where  Uj  is  given  by  (3.3). 

Step  7  If  U  —  0,  stop.  Otherwise,  remove  a  set  U  from  U  and  go  to  step  2. 

This  procedure  can  be  altered  in  the  usual  way  to  produce  bounds  (See  Section  2.6). 

A  few  comments  need  to  be  made  about  STDA-II.  We  first  look  at  Step  3.  When 
the  weights  of  the  arcs  are  raised  to  the  levels  j3,  the  tree  T(a)  will  not  in  general 
be  a  minimum  spanning  tree.  However,  if  its  weight  Wp{a)  does  not  exceed  d, 
then  by  definition  neither  will  W(0),  the  weight  of  a  minimum  spanning  tree  at  /3. 
Thus  in  such  a  situation  we  may  declare  U  feasible  without  actually  making  a  MST 
evaluation  at  0. 

If  W^(a)  >  d,  we  cannot  conclude  that  0  is  an  infeasible  point.  But  we  do  know 
that  there  is  a  point  with  lj  =  0j  for  j  0  T(a),  <  0j  for  at  least  one  j  £  T(a), 
at  which  T(a)  has  feasible  length.  (The  point  with  lj  =  otj  for  all  j  £  T(a)  is 
clearly  such  a  point  since  necessarily  W{ot)  <  d.)  Consistent  with  our  usual  goal  of 
keeping  the  list  of  undetermined  sets  small,  we  select  from  among  the  (potentially) 
many  such  points  l  one  which  minimizes  the  number  of  new  undetermined  intervals 
produced.  In  view  of  (3.3),  Uj  =  0  only  if  lj  =  0j  and  it  is  thus  reasonable  to  choose 
a  point  /  for  which  lj  =  0j  for  as  many  j  £  T(a)  as  possible.  Thus  we  have  a  0  —  1 
knapsack  problem  to  be  solved.  The  available  space  in  the  knapsack  is  d  —  W(a). 
Each  arc  j  £  T(a)  takes  up  space  equal  to  Wj{0j )  —  Wj(atj)  (here  putting  arc  j  into 
the  knapsack  means  setting  lj  =  0j)  and  all  arcs  have  equal  value,  say  1.  Steps  4 
and  5  apply  a  standard  knapsack  heuristic,  which  is  optimal  in  this  case  as  shown 
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below.  When  the  first  arc  is  found  that  “does  not  fit,”  we  increase  its  weight  as 
much  as  possible  from  the  level  ay  (i.e.,  “put  in”  as  much  of  it  as  possible)  and  leave 
the  remaining  tree  arcs  at  the  level  a. 

The  following  proposition  shows  that  the  steps  4  and  5  in  procedure  STDA-II 
result  in  a  maximum  number  of  arcs  in  T(a)  with  /y  =  /?y,  thus  producing  a  minimum 
number  of  new  undetermined  sets. 

Proposition  3.2  Steps  4  and  5  in  procedure  STDA-II  yield  a  minimum  number  of 
new  undetermined  intervals  in  each  iteration. 

Proof:  We  establish  this  result  by  showing  that  the  heuristic  knapsack  solution 
produced  by  these  steps  is  optimal.  Note  that  for  this  problem  all  items  considered 
for  inclusion  in  the  knapsack  (the  arcs  in  T(a))  have  equal  value.  Thus,  the  standard 
heuristic  merely  takes  items  in  increasing  order  of  u>j(/?y)  —  wJ(aJ). 

Consider  an  optimal  solution  with  value  n°.  If  there  is  any  arc  j  6  T(a)  not 
included  whose  weight  w3(0j)  —  Wy(ay)  is  less  than  that  of  any  arc  in  the  knapsack, 
these  arcs  can  clearly  be  exchanged  to  yield  a  solution  that  is  still  feasible  and  which 
also  has  value  n°  That  is,  this  produces  an  alternate  optimal  solution.  We  may 
repeat  this  procedure  until  all  arcs  in  the  knapsack  have  weights  at  least  as  small  as 
those  of  arcs  not  included,  again  yielding  an  optimal  solution.  But  this  solution  is 
precisely  the  solution  gotten  by  applying  the  above  heuristic.  This  proves  the  result. 
□ 

Finally,  since  only  those  arcs  j  E  T(a)  for  which  ay  <  /3y  are  considered  in  Steps 
4  and  5,  we  would  like  T(a)  to  contain  as  many  arcs  as  possible  for  which  ay  = 

(or  even  ay  +  1  =  /?y).  Thus  we  preprocess  the  arcs  prior  to  finding  T(a ): 

{iny(ay)  -  2e  if  ay  =  0, 

UJy(ay)  -  e  if  ay  +  1  =  /?y  , 

tny(ay)  otherwise 

where  t  is  suitably  small.  In  this  way  as  many  of  these  preferred  arcs  as  possible  are 
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included.  After  finding  T(a),  the  true  arc  weights  are  restored  so  that  tree  weights 
may  be  accurately  computed. 

STDA-II  is  an  intuitively  pleasing  algorithm  since  it  uses  easily  derived  points 
and  attempts  to  keep  the  list  of  undetermined  sets  small.  It  has  shown  promise  for  I 

fast  bound  production.  It  will  be  tested  on  several  examples  in  Section  3.4  along 
with  the  various  forms  of  STDA-I  and  ST-HYBRID. 

3.2.3  Hybrid  Spanning  Tree  Decomposition  Algorithm  1 

Thus  far  we  have  addressed  Problem  ST1  using  GSD  procedures.  We  describe 
here  ST-HYBRID,  an  algorithm  that  in  each  iteration  decomposes  an  undetermined 
set  U  =  [or,  0\  into  one  feasible  set  F,  one  infeasible  set  /  and  up  to  2a  —  1  new  I 

undetermined  sets.  Its  name  stems  from  the  fact  that  it  is  essentially  a  combina¬ 
tion  of  STDA-I  feasible  decomposition  and  STDA-I  infeasible  decomposition.  The 
following  description  summarizes  this  procedure. 

Procedure  ST-HYBRID 
Step  1  Start  with  U  =  fl.  Set  U  =  0,  Pf  =  0. 

I 

Step  2  Let  a  and  0  denote,  respectively,  the  lower  and  upper  endpoints  of  U.  Find 
T(0),  a  MST  at  0,  with  weight  W(0).  If  W(0)  <  d  then  U  is  feasible.  Set 
Pf  =  Pf  +  P{U}  and  go  to  step  7. 

I 

Step  3  Find  T(a),  a  MST  at  a,  with  weight  W’(a).  If  W{a)  >  d  then  U  is  infeasi¬ 
ble.  Go  to  step  7. 

Step  4  Derive  a  feasible  point  W  based  on  the  MST  T{a )  as  in  STDA-I  feasible 
decomposition. 

Step  5  Derive  an  infeasible  point  V  based  on  the  MST  T{0)  as  in  STDA-I  infeasible 
decomposition. 
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Step  6  Set  PF  =  PF  +  P{F},  where  F  is  defined  by  (3.2).  Decompose  the  remain¬ 
der  based  on  lF  and  V  as  described  in  the  proof  of  Proposition  2.3.  File  new 
undetermined  intervals  in  U. 

Step  7  If  14  =  0,  stop.  Otherwise  remove  a  set  U  from  U  and  go  to  step  2. 

This  procedure  can  be  altered  as  in  Section  2.6  to  produce  bounds. 

Section  2.5  showed  that  routines  like  ST-HYBRID  can  file  unnecessarily  large 
numbers  of  new  undetermined  sets  (compared  with  strict  GSD  routines).  It  was  also 
pointed  out,  though,  that  there  are  many  factors  that  jointly  determine  the  ultimate 
performance  characteristics  of  a  decomposition  routine.  In  Section  3.4  ST-HYBRID 
will  be  compared  with  the  other  routines  in  Section  3.2  on  some  sample  problems. 

We  close  this  discussion  of  Problem  ST1  with  the  following  GSD  modification 
that  simultaneously  determines  the  entire  distribution  of  W(X). 

3.2.4  Computing  the  Entire  Distribution  of  the  Minimum 
Spanning  Tree  Weight 

The  previous  sections  have  given  decomposition  algorithms  for  computing  a 
single  ordinate  of  the  cumulative  distribution  function  of  the  weight  of  a  MST. 
Calculating  P{  W (X)  <  d}  is  sufficient  for  situations  in  which  d  represents  an  ex¬ 
isting  constraint  on  MST  weight  (e.g.,  a  capital  or  material  constraint).  But  when 
that  is  not  the  case,  it  is  instructive  to  explore  the  full  range  of  possible  MST 
weights,  calculating  associated  probabilities  and  computing  the  mean  and  higher 
moments  of  W(X).  In  that  case,  we  could  execute  STDA-1  and/or  STDA-II  repeat¬ 
edly  for  all  possible  MST  weights  between  W(l,...,l)  and  W(»i, . . . ,  na).  This 
would  clearly  involve  a  great  deal  of  computational  effort,  particularly  when  there 
are  a  large  number  of  possible  realizations  of  W(X).  In  this  section  we  present  a 
GSD  modification  that  simultaneously  determines  the  entire  distribution  for  W(X) 
in  a  single  decomposition  procedure.  Procedure  STDIST  below  computes  the  pmf 
q(d)  =  P{W(X)  =  d).  Our  objective  requires  a  modification  to  the  definition  of 
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a  feasible  set.  For  each  undetermined  interval,  the  feasible  states  will  have  MST 
weights  equal  to  that  at  the  state  a.  (Clearly  all  variables  are  lower  feasible.)  We 
begin  by  finding  T(a),  a  MST  with  arc  weights  corresponding  to  the  point  a.  Prop¬ 
erties  3  and  4  yield  the  feasible  interval  defined  by  the  point 


r  _  f  <*j  3  €  T(a) 

J  ~  i  0}  if  j*T(a)  ■ 


(3.5) 


Up  to  a  undetermined  sets  also  .  esult  from  this  partition,  using  (3.3)  in  Section  3.2.2. 
This  procedure  is  summarized  below.  The  description  assumes  integer  arc  weights. 
Modifications  would  need  to  be  made  only  in  the  data  structures  for  the  q(d)'s  and 
the  manner  in  which  they  are  accumulated  if  arc  weights  are  non- integral. 


Procedure  STDIST 


Step  1  Find  W(l,...,l),  W(nj, . . .  ,na).  For  d  =  W(l, . . . ,  1), . . . ,  W(ni, . . .  ,n„), 
set  qi(d)  =  0.  Set  U  =  0  and  U  =  Q. 

Step  2  Let  a  and  /?  denote,  respectively,  the  lower  and  upper  endpoints  of  U.  Find 
T(a),  a  MST  at  a,  with  weight  W(a). 

Step  3  Let  l  be  as  defined  in  (3.5). 

Step  4  Put  <7f(W(a))  =  qi(W(a))  +  P{F},  where  F  is  defined  by  (3.2). 

Step  5  For  j  £  A:  If  lj  ^  fy,  file  Uj  in  U,  where  Uj  is  defined  by  (3.3). 

Step  6  If  U  —  0,  stop.  Otherwise,  remove  a  set  U  from  U  and  go  to  step  2. 

Note  that,  in  view  of  (3.5),  at  most  n  —  1  undetermined  sets  will  be  filed  in 
step  5.  Fewer  than  that  number  are  filed  if  a;  =  for  some  j  £  T(a).  For  that 
reason,  if  there  are  alternate  optimal  trees  at  a  we  would  like  to  choose  the  one  that 
contains  as  many  tree  arcs  as  possible  with  a,  =  0,.  In  order  to  accomplish  this, 
we  preprocess  the  arcs  in  a  manner  similar  to  that  used  in  STDA-II.  Specifically,  if 
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aj  =  /3j  we  replace  the  arc  weight  Wj(atj)  by  wj(aj)  —  e,  where  e  is  suitably  small. 

After  T(a )  is  found,  the  original  arc  weights  are  restored  so  that  VF(a)  may  be 
accurately  computed. 

As  usual,  STDIST  may  be  terminated  prior  to  complete  partitioning  of  fl  to  • 

produce  bounds.  To  produce  bounds  on  F(d)  =  /,{W(X)  <  d}  for  all  d ,  we 
substitute  step  6  for  step  6,  and  add  steps  7  and  8. 

Step  6  If  U  =  0,  stop.  If  a  stopping  condition  is  met,  go  to  step  7.  Otherwise,  9 

remove  a  set  U  from  U  and  go  to  step  2. 

Step  7  For  d  =  W(l, . . . ,  1), . . . ,  W(nu . . .  ,na):  let  F,(d)  =  Fn{d)  -  £?=w(i . 

Step  8  For  all  U  G  U:  Let  a  and  0  denote,  respectively,  the  lower  and  upper  ^ 

endpoints  of  U.  Find  VF(a)  and  W(0). 

For  W(a)  <d<  W(/3):  Set  Fu(d)  =  Fu{d)  +  P{U). 

For  W(/3)  <  d  <  W(n,\, . . .  ,na)  :  Set  Fu(d)  =  Fu{d)  +  P{U}  and  Ft{d)  =  9 

Fl{d)^P{U). 

3.3  Computing  Criticality  Indices  for  the 

Stochastic  MST  Problem  % 

We  now  describe  a  procedure  for  computing  PF{e)  =  P{ arc  e  is  on  a  MST},  the 
criticality  index  for  arc  e  as  given  in  Definition  1.1.  We  first  show  that  this  problem 

ft 

is  #P-hard. 

Proposition  3.3  The  evaluation  of  PF(e)  =  P{arc  e  is  on  a  MST}  is  #P~hard. 

Proof:  Let  G  =  (A/*,  A)  be  a  graph  whose  arcs  have  random  lengths  with  possible  I 

values  in  {0, 1}.  Let  s,t  G  A/”  and  let  X  be  the  vector  of  arc  weights,  which  corre¬ 
spond  here  to  lengths.  Let  L(X)  designate  the  (random)  length  of  a  shortest  s  —  t 
path  in  G  when  arc  lengths  are  given  by  X.  The  evaluation  of  P{L(\)  >  1}  is 
#P-hard  [3].  1 
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Now  create  an  arc  e  =  (s,  t)  with  fixed  length  1  and  let  G'  =  (Af,  A  +  e).  Then  * 

evaluating  P{arc  e  is  on  a  MST  in  G'}  is  equivalent  to  evaluating  P{L(\)  >  1}.  To 
see  this,  we  need  only  establish  that  for  a  given  realization  x  of  arc  lengths,  e  is  on 
a  MST  in  G'  if  and  only  if  L(x )  >  1.  This  is  proved  most  easily  by  contrapositive.  ^ 

Suppose  e  can  be  on  no  MST  in  G'.  Consider  a  MST  T,  and  let  C(e)  be  the 
unique  cycle  that  e  forms  with  T.  Every  arc  in  C(e)  —  e  must  have  zero  length,  for 
otherwise  e  could  be  pivoted  in  to  replace  an  arc  in  C(e)  —  e  with  length  1,  yielding 
an  MST  containing  e.  But  C(e)  —  e  necessarily  is  an  s  —  t  path  with  zero  length.  * 

Then  L(x)  =  0. 

Now  suppose  L(x)  —  0.  Consider  P,  an  s  —  t  path  in  G  of  length  0,  and  an 
arbitrary  MST  V  in  G' .  If  P'  =  P\T'  is  nonempty,  then  we  may  pivot  each  arc  ^ 

in  P'  into  T'  one  at  a  time,  arriving  at  a  new  MST  T  containing  P.  Since  T  is 
acyclic,  e  £  T.  Furthermore,  since  necessarily  C(e)  —  e  =  P,  e  cannot  be  pivoted  in 
to  achieve  a  MST  containing  e.  Then  no  MST  contains  e. 

Therefore,  * 

PF(e)  =  P{arc  e  is  on  a  MST}  =  P{L(X)  >  1} 

and  the  evaluation  of  PF(e )  is  #P-hard.  □  ) 

In  finding  PF(e),  the  feasible  region  consists  of  all  points  in  fl  at  which  arc  e 
is  on  a  MST.  Before  presenting  any  details  of  the  GSD  solution  of  this  problem, 
though,  we  need  to  investigate  the  nature  of  feasibility  in  this  setting. 

I 

Proposition  3.4  Let  /eft,  and  let  the  weight  of  arc  j  6  A  be  w3(J.j).  Suppose  that 
arc  e  is  in  a  MST  at  this  level.  Then  e  is  in  a  MST  for  all  states  in  the  set 

F  =  {leP  :  1  <  le  <  Je  (3.6)  > 

h  <  lj  <  n_,  for  j  /  c}. 

Proof:  All  points  in  F  can  be  arrived  at  from  /  by  a  sequence  of  steps,  each  involving 

* 

a  single  arc.  Possible  steps  involve  either  decreasing  the  weight  of  e  or  increasing 
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that  of  an  arc  j  ^  e  (either  j  G  T  or  j  £  T ).  Thus  we  need  only  show  that  each 
such  step  preserves  the  requirement  that  e  is  on  a  MST.  We  consider  each  of  three 
types  of  steps  separately. 

i)  Increasing  the  weight  of  arc  j  ^  e,j  G  T.  By  Property  4  (Section  3.1),  this 
step  either  does  not  affect  the  MST  T  or  causes  the  MST  to  change  to  T  —  j  +  i, 
where  i  £  T.  In  either  case,  since  j  ^  e,  e  can  stay  on  a  MST. 

ii)  Increasing  the  weight  of  arc  j  ^  e,  j  ^  T.  By  Property  3,  this  step  cannot 
force  a  change  in  T.  Thus  e  is  still  on  a  MST. 

iii)  Decreasing  the  weight  of  arc  e  G  T.  By  Property  1,  T  is  still  a  MST,  and 
eeT. 

Therefore,  if  /  is  feasible  then  F  is  a  feasible  set.  □ 

Proposition  3.5  Let  l  G  fl,  and  let  the  weight  of  arc  j  G  A  be  Wj(lj).  Suppose  that 
arc  e  is  in  no  MST  at  this  level.  Then  e  is  in  no  MST  for  any  point  in  the  set 

/  =  {/€//:  l  <lc<ne  (3.7) 

1  <  h  <  h  forj^e}. 

Proof:  All  points  in  /  can  be  arrived  at  from  /  by  a  sequence  of  steps,  each  involving 
a  single  arc.  Possible  steps  involve  either  increasing  the  weight  of  arc  e  or  decreasing 
the  weight  of  an  arc  j  ^  e  (either  j  G  T  or  j  T).  Thus  we  need  only  show  that 
each  such  step  preserves  the  requirement  that  e  can  be  on  no  MST.  Note  that  this 
property  is  equivalent  to  requiring  that  the  weight  of  arc  e  (strictly)  exceeds  the 
weights  of  all  arcs  in  C(e)  —  e,  where  C(e)  is  the  unique  cycle  that  e  forms  with  the 
MST  T.  We  consider  each  of  three  types  of  steps  separately. 

i)  Decreasing  the  weight  of  arc  j  ^  e,  j  G  T.  If  j  G  (7(e)  -  e,  then  the  weight  of 
the  most  heavily  weighted  arc  in  (7(e)  —  e  either  stays  the  same  or  is  decreased.  If 
j  &  (7(e)  —  e,  then  (7(e)  —  e  is  unchanged.  In  all  of  these  cases,  the  weight  of  arc 
e  still  exceeds  that  of  the  most  heavily  weighted  arc  in  C(e)  —  e,  so  that  e  is  in  no 
MST. 
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ii)  Decreasing  the  weight  of  arc  j  ^  e,  j  0  T.  By  Property  2,  this  step  either 
does  not  impact  T,  in  which  case  e  clearly  is  still  in  no  MST,  or  causes  an  arc  i  G  T 
to  be  displaced  from  T  by  arc  j.  If  i  ^  (7(e)  —  e,  the  cycle  that  e  forms  with  the 
new  MST  is  the  same  as  it  was  before,  namely  (7(e)  —  e,  so  e  still  cannot  pivot  in. 
If  i  €  (7(e)  —  e,  then  the  cycle  that  e  forms  with  the  new  MST  71  —  i  +  j  is  different 
from  (7(e)  —  e.  Specifically,  the  new  cycle  includes  the  remnants  of  the  old  cycle, 
(7(e)\{e,  *},  plus  the  remnants  of  the  cycle  j  formed  with  T,  C(j)  —  i.  The  weight 
of  e  still  exceeds  that  of  all  arcs  in  C(e)\{e,i}.  Furthermore,  it  must  be  that  arc  i 
was  the  most  heavily  weighted  arc  in  C(j)  —  j  (since  i  was  pivoted  out).  Since  we 
know  the  weight  of  e  exceeds  that  of  arc  i  (it  couldn’t  displace  it),  it  also  exceeds 
the  weights  of  all  arcs  in  C(j).  So  e  still  cannot  displace  any  arc  in  its  cycle  with 
the  new  MST,  therefore  it  can  be  on  no  MST. 

iii)  Increasing  the  weight  of  arc  e  T.  Clearly,  the  weight  of  arc  e  still  exceeds 
that  of  all  arcs  in  C(e)  —  e,  so  that  e  can  be  on  no  MST. 

Therefore,  if  /  is  infeasible  then  /  is  an  infeasible  set.  □ 

An  immediate  consequence  of  these  propositions  is  that  the  random  variable  Xt 
is  lower  feasible  while  the  remaining  variables  are  upper  feasible.  Thus  for  U  =  [a,  /?] 
the  extreme  feasible  candidate  is  the  point  (ae j/?)  =  (/?i, . . . ,  /3e_i ,  ae,  /?e+i, . . . ,  /?„) 
and  the  extreme  infeasible  candidate  is  ( fie\a )  =  (al5 . . . ,  ac_i,  ^e,ae+i, . . . ,  aa).  As 
with  other  GSD  applications,  the  above  propositions  also  imply  that  we  may  pre¬ 
screen  an  undetermined  set  U  —  [o,  /3]  by  evaluating  the  extremes.  If  (ae|/J)  is 
infeasible  then  U  is  infeasible,  and  if  (/?e|<*)  is  feasible  then  U  is  feasible. 

We  shall  now  describe  procedure  STCRIT.  As  in  STDIST,  limiting  indices  lj 
will  not  be  used.  Thus  decomposition  from  the  feasible  extreme,  in  which  we  push 
in  from  extreme  feasible  point  to  derive  a  feasible  point  /,  will  yield  one  feasible  set 

F=  {/  €U:  ae<le<Ie  (3.8) 

ii<U<Pi  fori^e}, 


I 
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the  undetermined  set 


Ue  =  {leU:  l+l  <l'<  pe  (3.9) 

Qi  <  h  <  Pi  for  i  ±  e}  (  =  0  if  le  —  /?«.), 

and  for  j  ^  e  the  undetermined  sets 

Uj  =  {/  G  U  :  <  /,  <  Pi  for  i  <  j,  i  ^  e 

1  (3.10) 

a,  <  f«  <  Pi  for  i  >  j,  i  /  e 
(*'  <le  <  !'  }  (  =  0  if  lj  =  ocj). 

Decomposition  from  the  infeasible  extreme,  in  which  we  push  in  from  the  extreme 
infeasible  point  to  derive  an  infeasible  point  /,  will  yield  one  infeasible  set 

I  =  {l€V:  ae<l'<  l  (3.11) 

a,  <  l,  <  h  for  i  /  e), 

the  undetermined  set 

Ue=  {l  £U:  ae<le<l'~  1  (3.12) 

a,  <  l,  <  Pi  for  i  ^  e]  (  --  0  if  lF  =  a,). 

and  for  j  ^  e  the  undetermined  sets 

=  {l  E  U  :  a,  <  l,  <  l,  for  i  <  j,  i  /  e 

h  +  l<h<  Pi  (3.13) 

a,  <  /,  <  p,  for  i  >  j ,  i  e 
l<lr<Pr  }  (  =  0  if 


As  in  STDA-II,  the  prescreening  phase  and  the  phase  in  which  we  decide  whether  to 
decompose  feasibly  or  infeasibly  are  somewhat,  free-flowing.  A  number  of  possible 
feasible  end  infeasible  points  are  investigated  first,  in  hope  of  deriving  a  particularly 


effective  and/or  easily  obtained  cutoff  state.  If  that  fails,  we  may  choose  between 
default  feasible  and  infeasible  decompositions.  Throughout,  Properties  1  through  4 
from  Section  3.1  will  be  used.  Steps  11  through  14  may  be  used  if  bounds  on  PF{e) 
are  desired.  We  present  this  procedure  in  narrative  form. 

Procedure  STCRIT(e) 

Step  1  Start  with  (J  =  fi.  Set  U  =  0  and  Pf  (e)  =  0. 

Step  2  Let  a  and  0  denote,  respectively,  the  lower  and  upper  endpoints  of  U.  Find 
T(0e\a),  a  MST  at  the  extreme  infeasible  candidate.  If  e  €  T((3e\a ),  go  to 
step  3.  Otherwise,  find  ij,  the  most  heavily  weighted  arc  on  the  cycle  that  arc 
e  forms  with  T(0e |a).  If  arc  e  and  arc  z'i  have  equal  weight,  pivot  arc  e  into 

T(0e\a). 

Step  3  If  e  ^  T(0e\a),  go  to  step  4.  Otherwise,  (/?e|a)  is  a  feasible  point  so  that  U 
is  feasible.  Go  to  step  1 1 . 

Step  4  If  U7e(ae)  >  te, ,(<*,, )  then  T(0e\a)  is  a  MST  at  levels  a  (call  it  T(a))\  go  to 
step  5.  Otherwise,  decreasing  the  weight  of  arc  e  will  cause  it  to  displace  arc 
i i.  Let  le  =  max {7  :  ae  <  7  <  0t\  wc{i)  <  lOq(orq)}  and  /”  =  a:  for  j  ^  e. 
Perform  feasible  decomposition  and  go  to  step  11. 

Step  5  Find  T(at\f)),  a  MST  at  the  extreme  feasible  candidate.  If  e  6  T(ac\0),  go 
to  step  6.  Otherwise,  find  f2,  the  most  heavily  weighted  arc  on  the  cycle  that 
e  forms  with  T(ae\f3).  If  arc  e  and  arc  i2  have  equal  weight,  pivot  arc  e  into 
T(ae \0)  and  go  to  step  6.  Otherwise,  e  is  on  no  MST  at  (af|/?),  which  is  thus 
an  infeasible  point.  Since  U  is  infeasible,  go  to  step  11. 

Step  6  We  have  r  £  T(ae\(3).  We  now  seek  to  force  e  out  of  the  MST  by  changing 
the  levels  of  only  a  few  arcs,  yielding  an  infeasible  point.  One  preprocessing 
pass  through  the  arcs  not  in  T(ae\0)  (similar  to  the  procedure  SHORT  in 
Section  3.2.1,  but  slightly  more  involved)  finds  the  following: 
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•  jx  =  the  least  weighted  arc  off  T(ae\ft)  whose  cycle  contains  arc  e,  which 
could  displace  arc  e  if  the  weight  of  e  is  increased  (in  step  7). 

•  j2  =  an  arc  not  in  T(ae\/3)  which  has  e  (at  level  /3)  as  its  most  heavily 
weighted  cycle  arc  and  for  which  Wj2{a32)  <  we{fir).  If  the  weight  of  e 
has  been  increased  to  tce(/?e)  (in  step  7),  then  j2  can  be  lowered  (in  step 
8)  to  force  e  out  of  the  tree.  If  no  such  arc  exists,  let  j2  =  0. 

•  j3  =  an  arc  not  in  T(ae\0 )  whose  cycle  with  T(ae  |/3)  contains  e  (not  as 
the  most  heavily  weighted  cycle  arc),  for  which  Wj{a3)  <  wt ( /3e )  for  all 
j  6  C{j3)  —  e  (for  use  in  step  9).  If  no  such  arc  is  found,  set  j3  —  0. 

Note:  Even  though  all  these  arcs  might  not  be  needed,  we  may  efficiently 
determine  all  three  in  one  O(an)  procedure,  essentially  the  same  amount  of 
work  needed  to  determine  any  one  of  these  numbers. 

Step  7  If  we(/3c)  <  wn(0n),  then  increasing  the  weight  of  arc  e  to  we(/3e)  does  not 
force  it  out  of  the  current  MST;  go  to  step  8.  Otherwise,  this  increase  causes 
arc  e  to  be  displaced  by  j\.  Set  le  =  min{7  :  ae  <  7  <  /3e  :  we(i)  >  u>j, (/?,,)} 
and  lj  =  (3j  for  j  ^  e.  Perform  infeasible  decomposition  and  go  to  step  1 1 . 

Step  8  T(ae\0)  is  now  a  MST  at  levels  /?;  call  it  T(f3).  If  j2  =  g°  to  step  9. 

Otherwise,  decreasing  the  weight  of  j2  will  cause  it  to  displace  e  from  T(j3). 

Set  /e  =  &,  ll2  -  max {7  :  an  <7  <  0j2;wJ2{~i)  <  u>e(&)}  and  l3  =  for  all 
other  arcs  j.  Perform  infeasible  decomposition  and  go  to  step  11. 

Step  9  If  j3  =  0,  go  to  step  10.  Otherwise,  for  j  £  C{j3)  —  j3  —  e,  set  l3  =  max{7  : 
a}  <  7  <  <  u?e(/3e)}  and  =  f3e.  Then  e  is  the  most  heavily- 

weighted  arc  on  the  cycle.  Set  ln  =  max{7  :  an  <  7  <  /J73;  ^(7)  <  wr ( fie ) } 

to  displace  e.  Set  l3  =  for  all  arcs  j  whose  weights  are  not  yet  set.  Perform 
infeasible  decomposition  and  go  to  step  1 1 . 

Step  10  We  are  left  to  choose  from  the  default  decompositions.  Do  one  of  the 
following: 
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•  Derive  a  feasible  point  from  the  known  feasible  point  0  by  setting  l3  =  a3 
for  j  G  T(0),  le  =  0e,  and  lj  =  min{7  :  otj  <  7  <  0j\j  does  not  displace  e 
from  T(0)}  for  j  ^  T(0).  Perform  feasible  decomposition  and  go  to  step 
11. 

•  Derive  an  infeasible  point  from  the  known  infeasible  point  a  (which  we 
know  to  be  infeasible  from  step  4  above)  by  setting  lj  =  03  for  j  T(a), 
le  =  Qfe,  and  lj  —  max{7  :  oij  <  7  <  03;j  stays  in  T(a)  when  arcs 
j  T(a)  are  at  levels  ct,  }  for  j  G  T(a).  Perform  infeasible  decomposition 
and  go  to  step  1 1 . 

Step  11  If  one  or  more  feasible  intervals  have  been  generated,  add  their  probability 
to  PtF (e).  If  U  =  0,  stop.  If  a  stopping  condition  is  met,  go  to  step  12. 
Otherwise,  remove  a  set  U  from  U  and  go  to  step  2. 

Step  12  Set  PF(e)  =  PtF(e) 

Step  13  For  all  U  G  U\  Let  a  and  0  denote,  respectively,  the  lower  and  upper 
limiting  points  of  U .  Find  T(ae \0).  If  e  G  T(ac\0)  or  if  e  can  be  pivoted  into 
T(ae\0)  yielding  an  alternate  MST,  set  PF(e)  =  PF(e)  +  P{U}. 

Step  14  Return  bounds  P/(e)  <  PF(e)  <  P„(e). 

A  few  comments  are  in  order  concerning  STCRIT.  In  steps  3  through  9  several 
attempts  are  made  to  find  feasible  or  infeasible  points  that  are  easy  to  find  and  that 
will  yield  few  undetermined  sets.  If  U  is  found  to  be  all  feasible  in  step  3,  clearly  no 
undetermined  sets  are  produced.  If  reducing  arc  e  in  step  4  yields  a  feasible  point, 
then  one  undetermined  set  will  be  produced.  Finding  U  to  be  infeasible  in  step  5 
produces  no  undetermined  sets.  If  raising  arc  e  from  (ae|/3)  causes  it  to  pivot  out  of 
the  tree  in  step  7,  then  one  undetermined  set  is  produced.  If  arc  j2  can  be  lowered 
so  as  to  displace  arc  e  from  T(0)  in  step  8,  two  undetermined  sets  are  produced. 
Finally,  if  step  9  succeeds  in  identifying  an  infeasible  set  we  get  undetermined  sets 
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for  arc  e,  for  j3  and  for  each  of  the  arcs  on  the  cycle  formed  by  j3  with  T(jd).  If  we 
are  left  with  the  default  decompositions  of  step  10,  we  get  up  to  n  —  1  undetermined 
sets  (if  we  decompose  feasibly)  or  up  to  a  —  n  +  1  (if  we  decompose  infeasibly).  In 
that  case,  the  decomposition  yielding  fewest  new  undetermined  sets  is  used. 

One  final  comment  about  STCRIT  is  that  it  can  be  used  with  little  modification 
in  a  way  similar  to  that  given  in  [l]  to  compute  the  P{arc  e  is  on  a  MST  and 
IU(X)  >  d}.  Combining  this  measure  with  the  output  of  STDA-I  or  STDA-II,  we 
can  compute  the  conditional  criticality  index 


P{arc  e  is  on  a  MST  |  W(X)  >  d}  = 


arc  e  is  on  a  MST  and  W(. 

1  -P{W(X)<d\ 


These  probabilities  are  useful  in  that  they  identify  the  arcs  that  are  likely  to  be 
MST  members  when  the  network  cannot  be  connected  feasibly.  Unfortunately,  it 
takes  a  bit  of  extra  work  to  produce  them.  We  first  run  STDA-I  or  STDA-II,  filing 
all  generated  infeasible  sets  in  a  list  1.  If  these  are  run  to  complete  decomposition, 
the  union  of  the  sets  in  I  is  precisely  the  infeasible  region  of  D.  These  infeasible 
sets  are  passed  as  input  to  STCRIT  (rather  than  starting  with  U  =  0  in  step  I,  we 
begin  with  U  =  T  and  instead  of  using  fl  as  the  first  undetermined  set  we  simply 
choose  U  from  U). 

In  conclusion,  STCRIT  is  an  algorithm  which  can  produce  marginal,  joint  or 
conditional  criticality  indices  for  any  arc  e  E  A.  These  measures  are  useful  in 
determining  which  arcs  are  influential  in  the  formation  of  minimum  spanning  trees. 
Section  3.4  below  contains  examples. 


3.4  Stochastic  MST  Examples 

In  this  section  we  demonstrate  the  performance  of  the  various  routines  developed 
in  the  present  chapter.  For  each  example,  we  execute  STDA-I  (Pure  Strategies  1 
and  2,  Mixed  Strategies  1  through  3),  STDA-II  (Pure  Strategy  1  only),  STDIST 
and  STCRIT.  Our  main  goals  are  to  compare  the  performance  of  several  versions 
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of  STDA-I  and  STDA-II,  and  to  show  examples  of  complete  decomposition  as  com¬ 
pared  to  partial  decomposition  to  produce  bounds.  All  codes  were  implemented 
in  FORTRAN  77  and  were  compiled  and  executed  on  a  SUN  SPARCstation  IPC. 
The  implementation  of  Kruskal’s  algorithm  in  Camerini,  Galbiati  and  Maffioli  [12] 
was  used  for  computing  a  minimum  spanning  tree.  In  all  applications,  the  list  of 
undetermiend  sets  was  kept  in  a  heap  sorted  by  probability  until  the  most  probable 
undetermined  set  had  probability  less  than  10~2°.  Arc  weight  distributions  repre¬ 
sent  what  we  feel  to  be  reasonable  situations.  Specifically,  a  high  probability  is 
assigned  to  the  lowest  weight  (the  status  quo)  and  low  probabilities  are  assigned  to 
the  higher  weights  (the  occasional  deviations  from  the  status  quo).  All  CPU  times 
quoted  are  in  seconds.  In  all  tables,  P{U}  represents  the  sum  of  the  probabilities 
of  undetermined  sets  contained  in  the  collection  U. 

Example  3.1  (Example  1.1  revisited)  The  network  in  Table  3.1  has  10  nodes  and 
21  arcs.  Each  arc  has  three  possible  weights.  The  state  space  has  more  than 
10  billion  states.  For  this  example,  the  range  of  possible  MST  weights  is  from 
W(l,...,l)  =  47to  W(3,...,3)  =  138.  We  use  STDA-I,  STDA-II  and  ST-HYBRID 
to  compute  P{W(X)  <  60}  and  F{W(X)  <  90},  these  tree  weights  being  relatively 
extreme  possible  values,  chosen  to  demontrate  the  difference  in  feasible  and  infeasible 
decompositions.  The  results  for  F(60)  are  in  Table  3.2,  for  F(90)  in  Table  3.3.  A 
partial  listing  of  the  cumulative  distibution  function  F(d)  (from  STDIST)  !s  given 
in  Table  3.4.  Full  decomposition  used  445.91  CPU  seconds  and  decomposed  634,405 
sets.  Partial  decomposition  decomposed  100,000  intervals  in  69.63  CPU  seconds. 
Criticality  indices  (from  STCRIT)  are  in  Table  3.5.  Arc  10,  which  stands  out  as 
requiring  more  computational  effort  than  other  arcs,  achieved  bounds  of  0.0001852 
and  0.0003542  in  6.81  seconds  after  decomposing  500  intervals. 

There  are  a  few  items  worthy  of  comment  in  the  GSD  solutions  of  Example  3.1. 
First,  we  note  the  incredible  efficiency  of  the  procedure  STCRIT.  Note  that  almost 
all  arcs  required  the  evaluation  of  fewer  than  200  undetermined  intervals.  Since  each 
such  evaluation  entails  two  MST  evaluations,  most  arcs  required  the  evaluation  of 
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Table  3.1:  Input  Distribution  for  Example  3.1. 


Arc 

Weight  (Probability) 

1  =  (1,2) 

2.0  (0.90) 

8.0  (0.08) 

12.0  (0.02) 

2  =  (1,3) 

10.0  (0.85) 

24.0  (0.12) 

35.0  (0.03) 

3  =  (1,4) 

6.0  (0.88) 

18.0  (0.10) 

24.0  (0.02) 

4  =  (2,5) 

17.0  (0.75) 

35.0  (0.20) 

50.0  (0.05) 

5  =  (2,6) 

3.0  (0.68) 

7.0  (0.25) 

10.0  (0.07) 

6  =  (3,2) 

12.0  (0.85) 

22.0  (0.11) 

30.0  (0.04) 

7  =  (3,7) 

10.0  (0.80) 

19.0  (0.14) 

24.0  (0.06) 

8  =  (3,8) 

4.0  (0.75) 

6.0  (0.17) 

10.0  (0.08) 

9  =  (4,3) 

3.0  (0.84) 

7.0  (0.13) 

12.0  (0.03) 

10  =  (4,9) 

21.0  (0.85) 

33.0  (0.09) 

45.0  (0.06) 

11  =  (5,7) 

8.0  (0.77) 

14.0  (0.13) 

18.0  (0.10) 

12  =  (5,10) 

15.0  (0.76) 

21.0  (0.22) 

25.0  (0.02) 

13  =  (6,3) 

5.0  (0.65) 

10.0  (0.23) 

12.0  (0.12) 

14  =  (6,5) 

18.0  (0.94) 

27.0  (0.05) 

36.0  (0.01) 

15  =  (6,7) 

25.0  (0.80) 

38.0  (0.12) 

50.0  (0.08) 

16  =  (7,8) 

20.0  (0.87) 

30.0  (0.09) 

40.0  (0.04) 

17  =  (7,10) 

4.0  (0.75) 

9.0  (0.14) 

15.0  (0.11) 

18  =  (8,4) 

9.0  (0.82) 

13.0  (0.15) 

20.0  (0.03) 

19  =  (8,9) 

15.0  (0.79) 

28.0  (0.15) 

45. .0  (0.06) 

20  =  (9,7) 

10.0  (0.95) 

17.0  (0.03) 

30.0  (0.02) 

21  =  (9,10) 

8.0  (0.85) 

14.0  (0.10) 

25.0  (0.05) 

about  400  states,  or  a  fraction  3.8  X  10-8  of  the  state  space.  Even  arc  10,  whose 
criticality  index  required  the  evaluation  of  many  more  states,  only  looked  at  the 
fraction  8.0  x  10-6  of  the  state  space. 

Secondly,  we  comment  on  the  bounds  produced  by  STDIST  for  Table  3.4.  Clearly, 
the  bounds  corresponding  to  small  MST  lengths  (small  values  of  d)  converge  more 
quickly  than  those  for  large  values  of  d.  This  is  due  to  the  manner  in  which  these 
probabilities  are  accumulated,  starting  at  the  feasible  extreme  and  moving  up. 

Example  3.2  We  consider  again  the  network  in  Example  3.1,  but  with  revised 
probability  structure.  Here,  each  arc  has  two  possible  weights,  as  shown  in  Table  3.6. 
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Table  3.2:  Computing  F( 60)  =  0.8267495828  for  Example  3.1. 


Routine 

Decor 

Sets 

Full 

nposition 

Time 

500  Sets  Decc 
Bounds 

imposec 

Time 

m 

PM 

STDA-I  PI 

790 

5.63 

.8197061589 

.8278633467 

3.67 

192 

.008157 

STDA-I  P2 

7886 

68.31 

.6655013736 

.9071519952 

4.67 

620 

.2417 

STDA-I  Ml 

648 

5.86 

.8231632665 

.8273711991 

4.61 

123 

.004208 

STDA-I  M2 

5582 

48.56 

.6614555582 

.8989402274 

4.97 

585 

.2375 

STDA-I  M3 

4822 

40.41 

.7109925528 

.8718820125 

5.02 

689 

.1609 

STDA-II  PI 

1421 

1.30 

.8188268475 

.8298433389 

0.67 

128 

.01102 

ST-HYBRID 

2660 

10.24 

.8065717649 

.8416536725 

3.01 

324 

.03508 

Table  3.3:  Computing  F(90)  =  0.9999950999  for  Example  3.1. 


Routine 

F 

Decom 

Sets 

nil 

sosition 

Time 

500  Sets  D 
Bounds 

ecompo 

Time 

sed 

\U\ 

P[U) 

STDA-I  PI 

247853 

1603.52 

.94844450 

1.0000000 

3.39 

516 

MMi 

STDA-I  P2 

205576 

1655.32 

.98419966 

.99999965 

5.09 

2506 

ESI 

STDA-I  Ml 

88857 

727.00 

.98871880 

.99999990 

4.89 

1315 

.01128 

STDA-I  M2 

171400 

1389.58 

.96960607 

.99999991 

4.53 

647 

STDA-I  M3 

195337 

1589.95 

.98216012 

.99999966 

5.54 

2329 

.01784 

STDA-II  PI 

127683 

95.18 

.99899010 

.99999999 

0.76 

686 

.001010 

ST-HYBRID 

506168 

2115.43 

.90042352 

.99999991 

2.37 

573 

.09958 

Table  3.4:  Partial  listing  of  the  cdf  F(d)  for  Example  3.1. 


Full  Decomposition 
d  |  F(i) 

100,000  Sets  Decomposed 
Lower  j  Upper 

0.0984162 

0.0984161 

0.0984161 

0.1450503 

0.1450503 

0.1450503 

0.1838571 

0.1838571 

0.1838571 

0.2594979 

0.2594979 

0.2594979 

0.2790302 

0.2790302 

0.2790302 

52.0 

0.3801882 

0.3801882 

0.3801882 

53.0 

0.4426021 

0.4426021 

0.4426021 

54.0 

0.5169705 

0.5169705 

0.5169705 

55.0 

0.5817478 

0.5817478 

0.5817478 

56.0 

0.6376404 

0.6376404 

0.6376404 

57.0 

0.6932358 

0.6931849 

0.6932502 

58.0 

0.7424916 

0.7423291 

0.7425250 

59.0 

0.7888656 

0.7885259 

0.7889317 

60.0 

0.8267496 

0.8260885 

0.8268742 

61.0 

0.8621487 

0.8609763 

0.8623779 

62.0 

0.8894412 

0.8873433 

0.8898221 

63.0 

0.9126231 

0.9093929 

0.9132196 

64.0 

0.9316627 

0.9270028 

0.9325470 

65.0 

0.9468859 

0.94022.°9 

0.9481127 

70.0 

0.9874525 

0.9684075 

0.9909956 

75.0 

0.9976452 

0.9736052 

0.9994587 

80.0 

0.9996402 

0.9782001 

0.9999404 

85.0 

0.9999545 

0.9834824 

0.9999924 

90.0 

0.9999951 

0.9870540 

0.9999992 

100.0 

0.9999999 

0.9945591 

0.9999999 
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Table  3.5:  Criticality  indices  for  Example  3.1 


Arc 

Criticality  Index 

#  Sets 

Time 

1 

0.9393247 

12 

2 

0.0432205 

11 

\ 

3 

0.5859462 

6 

4 

0.0490331 

3199 

18.94 

5 

0.8227428 

200 

1.18 

6 

0.0021765 

8 

0.047 

7 

0.8004465 

77 

0.46 

8 

0.9353420 

198 

1.17 

9 

0.9084087 

9 

0.053 

10 

0.0001994 

41939 

248.28 

11 

0.9067588 

174 

1.03 

12 

0.0833359 

9 

0.053 

13 

0.6990894 

16 

0.095 

14 

0.0099071 

169 

1.00 

15 

0.0000000 

1 

0.006 

16 

0.0001650 

12 

0.071 

17 

0.9007350 

4 

0.024 

18 

0.0764335 

10 

0.059 

19 

0.1628822 

13 

0.077 

20 

0.2318770 

10 

0.059 

21 

_ 

0.8654815 

8 

. ,  .  1 

0.047 

Table  3.6:  Input  Distribution  for  Example  3.2. 


Arc 

Weight  (Probability) 

1  =  (1,2) 

70.0  (0.95) 

94.0  (0.05) 

2  =  (1,3) 

25.0  (0.80) 

52.0  (0.20) 

3  =  (1,4) 

42.0  (0.70) 

61.0  (0.30) 

4  =  (2,5) 

26.0  (0.85) 

50.0  (0.15) 

5  =  (2,6) 

58.0  (0.70) 

95.0  (0.30) 

6  =  (3,2) 

15.0  (0.8®) 

43.0  (0.20) 

-a 

IJ 

'oo 
— ■ j 

65.0  (0.60) 

75.0  (0.40) 

oo 

II 

'So 

oo 

59.0  (0.70) 

98.0  (0.30) 

9  =  (4,3) 

21.0  (0.90) 

68.0  (0.10) 

10  =  (4,9) 

89.0  (0.85) 

96.0  (0.15) 

11  =  (5,7) 

32.0  (0.80) 

67.0  (0.20) 

12  =  (5,10) 

63.0  (0.75) 

99.0  (0.25) 

13  =  (6,3) 

66.0  (0.80) 

98.0  (0.20) 

14  =  (6,5) 

60.0  (0.90) 

88.0  (0.10) 

15  =  (6,7) 

32.0  (0.75) 

48.0  (0.25) 

16  =  (7,8) 

16.0  (0.7) 

42.0  (0.30) 

17  =  (7,10) 

56.0  (0.80) 

71.0  (0.20) 

18  =  (8,4) 

90.0  (0.95) 

96.0  (0.05) 

19  =  (8,9) 

17.0  (0.80) 

45.0  (0.20) 

20  =  (9,7) 

3.0  (0.60) 

15.0  (0.40) 

21  =  (9,10) 

50.0  (0.75) 

66.0  (0.25) 

For  this  example,  the  possible  MST  weights  range  from  220  to  444.  The  cardinality 
of  the  state  space  is  just  over  2  million.  In  Table  3.7,  we  use  STDA-I,  STDA-II  and 
ST-HYBRID  to  compute  F(280),  and  in  Table  3.8  we  give  results  for  F(400).  A 
partial  listing  of  the  cdf  F(d)  is  given  in  Table  3.9.  Full  decomposition  required  6.96 
CPU  seconds  to  decompose  9,743  intervals.  Table  3.10  contains  criticality  indices. 

Again,  we  see  the  efficiency  of  STCRIT,  which  used  just  539  interval  evaluations 
total  to  compute  all  21  criticality  indices.  Most  arcs  required  the  evaluation  of  fewer 
than  20  states,  a  fraction  9.5  x  10~6  of  the  state  space. 
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Table  3.7:  Computing  F(280)  =  0.8567535213  for  Example  3.2. 


Routine 

Full 

Sets 

Decomposition 

Time 

STDA-I  Pure  1 

740 

3.93 

STDA-I  Pure  2 

356 

2.39 

STDA-I  Mixed  1 

289 

2.46 

STDA-I  Mixed  2 

336 

3.03 

STDA-I  Mixed  3 

295 

2.63 

STDA-II  Pure  1 

658 

0.62 

ST-HYBRID 

573 

3.20 

Table  3.8:  Computing  F(40Q)  =  0.9999977607  for  Example  3.2. 


Routine 

Full 

Sets 

Decomposition 

Time 

STDA-I  Pure  1 

740 

3.93 

STDA-I  Pure  2 

78 

0.76 

STDA-I  Mixed  I 

176 

1.39 

STDA-I  Mixed  2 

143 

1.16 

STDA-I  Mixed  3 

109 

0.99 

STDA-II  Pure  1 

465 

0.46 

ST-HYBRID 

462 

2.20 
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Table  3.9:  Partial  listing  of  the  cdf  F(d)  for  Example  3.2. 


d 

F(d) 

d 

F(d) 

220.0 

0.0925345 

275.0 

0.8182486 

221.0 

0.1242605 

280.0 

0.8567536 

225.0 

0.1242605 

226.0 

0.1489364 

285.0 

0.8871976 

227.0 

0.1573967 

290.0 

0.9202321 

231.0 

0.1573967 

232.0 

0.2190863 

300.0 

0.9554863 

233.0 

0.2448637 

234.0 

0.2464500 

310.0 

0.9774270 

235.0 

0.2464500 

236.0 

0.2788370 

320.0 

0.9892763 

237.0 

0.3061347 

238.0 

0.3281373 

330.0 

0.9953135 

239.0 

0.3337775 

240.0 

0.3337775 

340.0 

0.9979105 

241.0 

0.3409746 

242.0 

0.3516675 

350.0 

0.9991380 

243.0 

0.3588058 

244.0 

0.3766160 

360.0 

0.9997070 

245.0 

0.3852992 

246.0 

0.4104818 

370.0 

0.9998960 

247.0 

0.4297512 

248.0 

0.4810531 

380.0 

0.9999668 

249.0 

0.5093677 

250.0 

0.5187621 

390.0 

0.9999884 

255.0 

0.5673748 

400.0 

0.9999978 

260.0 

0.6490759 

410.0 

0.9999997 

265.0 

0.7218336 

420.0 

0.9999999 

270.0 

0.7599369 

440.0 

1.0000000 

ft 


ft 
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Table  3.10:  Criticality  indices  for  Example  3.2. 


Arc 

Criticality  Index 

#  Sets 

Time 

1 

0.0000000 

1 

0.0066 

2 

0.8740000 

4 

0.0263 

3 

0.2260000 

5 

0.0329 

4 

1.0000000 

1 

0.0066 

5 

0.1400000 

3 

0.0197 

6 

1.0000000 

1 

0.0066 

7 

0.0003105 

33 

0.2168 

8 

0.0420000 

4 

0.0263 

9 

0.9000000 

2 

0.0131 

10 

0.0000000 

1 

0.0066 

11 

0.8000360 

8 

0.0526 

12 

0.0387825 

74 

0.4862 

13 

0.0001656 

46 

0.3022 

14 

0.0162000 

5 

0.0329 

15 

1.0000000 

1 

0.0066 

16 

0.7600000 

3 

0.0197 

17 

0.2000000 

3 

0.0197 

18 

0.0000000 

1 

0.0066 

19 

0.2400000 

3 

0.0197 

20 

1.0000000 

1 

0.0066 

21 

0.7625270 

339 

2.2272 

ft 


ft 


ft 


ft 


ft 


ft 


ft 
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3.5  Extension  of  MST  Results  to  Other  Stochas¬ 
tic  Matroid  Settings 

The  GSD  applications  given  in  this  chapter  were  successful  due  in  no  small 
measure  to  the  matroid  structure  of  the  MST  problem.  Note  that  spanning  trees 
are  bases  for  the  graphic  matroid  with  ground  set  A  and  independent  sets  being 
acyclic  subgraphs.  The  following  matroid  properties  were  instrumental  in  making 
these  applications  more  effective: 

•  Matroid  problems  can  be  solved  by  the  greedy  algorithm  (here,  Kruskal's 
method),  so  that  evaluating  a  state  point  (finding  an  optimal  matroid  base  for 
a  given  weighting  of  the  elements  of  the  ground  set)  can  be  done  efficiently. 

•  When  one  element  of  the  ground  set  changes  weight,  it  is  relatively  easy  to 
determine  whether  a  given  base  is  still  optimal.  Furthermore,  by  the  exchange¬ 
ability  property  of  matroids  it  is  easy  to  form  a  new  optimal  base  by  pivoting, 
assuming  some  means  exist  for  finding  the  replacement  element.  This  is  im¬ 
portant  because  the  alternative  is  to  solve  a  new  optimization  problem  every 
time  an  element  changes  weight. 

Since  all  matroid  problems  share  these  features  to  an  extent,  it  is  reasonable  to 
assume  that  there  are  other  stochastic  matroid  problems  that  will  admit  efficient 
solution  by  CISI).  The  second  point  is  especially  important.  In  our  MST  work,  every 
exchange  of  elements  was  accomplished  by  investigating  the  unique  cycle  (minimal 
dependent  set)  formed  by  an  arc  not  on  a  MST  (an  element  of  the  ground  set  not  in 
the  optimal  base)  with  the  MST  (the  elements  making  up  the  optimal  base).  These 
comparisons  resulted  in  efficient,  unequivocal  determination  of  the  elements  to  be 
exchanged.  All  matroid  problems  share  this  feature  to  some  degr<v.  In  tin-  gen 
eral  (non- matroid)  case,  it  is  not  always  easy  to  determine  whether  weight  changes 
cause  a  solution  to  cease  being  optimal,  or  to  remedy  the  situation  by  identifying  the 
elements  to  be  exchanged.  In  the  absence  of  some  means  for  making  these  determi- 
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natons,  every  change  in  element  weight  will  necessitate  finding  an  optimal  solution 
from  scratch. 

3.6  Conclusions 

The  GSD  routines  presented  in  this  chapter  represent  real  gains  in  our  ability  to 
characterize  the  behavior  of  minimum  spanning  trees  in  stochastic  networks.  They 
were  particularly  effective  due  to  the  attractive  properties  of  matroids.  We  close 
this  chapter  with  some  observations  about  the  computational  results  of  Section  3.4 
and  some  directions  for  further  study. 

It  seems  clear  from  Tables  3.2,  3.3,  3.7  and  3.8  that,  despite  its  sophistication, 
STDA-I  cannot  compete  with  STDA-II  in  terms  of  speed.  Based  on  our  experience, 
it  seems  to  be  the  case  with  GSD  applications  that  frequently  a  simpler,  clever  de¬ 
composition  scheme  performs  better  than  one  that  labors  to  cut  off  as  much  of  an 
interval  as  possible  in  each  iteration.  In  particular,  even  though  the  mixed  strategies 
of  STDA-I  improved  upon  the  STDA-I  pure  strategies,  they  do  not  compare  favor¬ 
ably  with  STDA-II  on  the  sample  problems.  It  is  possible  that  mixed  strategies  will 
work  well  with  other  GSD  applications.  It  also  seems  clear  that  the  hybrid  routine 
did  not  outperform  the  best  of  the  GSD  routines  either  in  terms  of  computational 
effort  or  in  terms  of  achieving  tight  bounds. 

Secondly,  it  is  clearly  beneficial  in  STDIST  to  accumulate  probabilities  for  all 
possible  MST  weights  at  the  same  time.  To  see  this,  note  that  in  Example  3.1  the 
fastest  routine  needed  95  seconds  to  compute  F{ 90).  STDIST  computed  the  entire 
distribution  in  an  average  of  4.85  seconds  per  possible  MST  weight. 

Finally,  even  among  MST  routines,  the  procedure  STCRIT  stands  out  as  a 
very  efficient  routine.  It  would  be  good,  then,  to  find  applications  for  this  effective 
procedure.  For  example,  a  reasonable  heuristic  solution  to  the  probl  m  of  ident  ifying 
the  most  probable  minimum  spanning  tree  is  suggested.  After  identifying  criticality 
indices  for  all  arcs,  we  build  a  maximally  acyclic  connected  graph  by  considering 


arcs  one  at  a  time  in  decreasing  order  of  their  criticality  index,  discarding  any  arc 
that  would  form  a  cycle.  Investigating  the  effectiveness  of  this  proposed  heuristic 
seems  worthy  of  future  research. 
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CHAPTER  4 


The  Stochastic  Minimum  Cost 
Network  Flow  Problem 


4.1  Introduction 

An  important  problem  in  the  area  of  network  optimization  is  that  of  satisfying 
demand  for  a  commodity  by  using  a  least-cost  subset  of  transportation  arcs,  with 
flow  constrained  on  each  arc  by  a  finite  upper  bound.  This  problem  is  known  as  the 
(capacitated)  minimum  cost  flow  (MCF)  problem. 

In  this  chapter  we  consider  flow  networks  in  which  neither  the  cost  of  using  an 
arc  nor  the  capacity  of  that  arc  is  known  with  certainty.  Rather,  these  values  are 
given  as  discrete  random  variables.  After  some  formal  preliminaries,  we  shall  present 
applications  of  the  General  State  Space  Decomposition  (GSD)  algorithm  of  Chapter 
2  designed  to  assist  in  characterizing  the  probabilistic  behavior  of  this  system. 

Consider  a  network  G  =  (AT,  A,s,t)  with  node  set  Af  =  {l,...,n}  and  arc  set 
A  =  { 1 , . . . ,  a}.  From  Af,  we  designate  nodes  s  and  t  as  the  unique  source  (supply  of 
the  commodity)  and  sink  (demand  for  the  commodity),  respectively.  Of  course,  this 
is  not  restrictive  since  the  existence  of  multiple  supply  and/or  demand  locations  can 
be  modeled  by  creating  a  fictitious  supersource  s  and  a  fictitious  supersink  t.  We 
assume  that  supply  at  s  equals  demand  at  t,  both  given  by  v  >  0.  We  assume  that 
nodes  do  not  restrict  flow  and  that  flow  passes  through  nodes  without  incurring  a 
cost.  Nodes  for  which  these  assumptions  do  not  hold  may  be  split  into  two  nodes 
connected  by  an  arc  with  the  desired  cost  and  capacity  properties. 
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Associated  with  each  arc  j  6  A  is  a  cost  X3  per  unit  of  flow  transmitted,  a 
discrete  random  variable  of  the  form  given  in  Chapter  2.  That  is,  the  cost  of 
using  arc  j  takes  on  finite  values  Cj(l)  <  •  •  •  <  Cj(mj)  with  respective  probabilities 
Pj(  1), . . .  Furthermore,  the  random  capacity  (upper  bound)  of  arc  j,  given 

by  Yj,  takes  on  nonnegative  finite  values  Uj(l)  <  •••  <  Uj(rtj)  with  probabilities 
^j(l), . . .  ,qj(rij).  We  represent  the  collection  of  cost  and  capacity  random  variables 
in  vector  form  as  (X|Y)  =  (Xj, . . . ,  Xa\Yi, . . . ,  Ya).  The  state  space  Q  of  X  contains 
n;=t  m  j  11?=  i  n3  points  of  the  form  (x|t/)  =  (c,(fc,), . . .  ,c0(A:a)|u1(Z1), . . . ,  ua(la)), 
where  kj  €  { 1 , . . . ,  m2  }  and  lj  E  { 1 ,  • .  • ,  n} }  for  all  j  £  A.  As  in  Chapter  2,  we  shall 
use  the  convention  of  writing  the  state  point  (x|y)  as  (fc|Z)  =  (jfcj, . . . ,  fca|/i, . . . ,  /„). 
For  a  fixed  (x|y)  6  fi,  let  C(x|y)  be  the  cost  of  a  MCF  (with  value  v). 

We  consider  first  the  following  problem  related  to  the  stochastic  MCF  system 
described  above. 

Problem  MCI  Let  d  >  0.  Assuming  the  components  of  (X|Y)  are  independent, 
evaluate  the  probability  P{C(X|Y)  <  d}  that  there  is  a  flow  of  v  units  whose 
cost  is  at  most  d. 

For  a  given  d,  a  realization  (k\l)  is  considered  to  be  a  feasible  state  if  C(k\l)  <  d. 
Clearly,  feasible  points  correspond  to  low  arc  costs  and  high  arc  capacities.  Thus 
the  cost  variables  Xj  are  lower  feasible  and  the  capacity  variables  Y3  are  upper 
feasible.  In  Section  4.2  we  compute  P{C(X|Y)  <  d]  and  give  a  GSD  modification 
that  simultaneously  determines  the  entire  probability  distribution  of  C(X|Y). 

In  Problem  MC2  below,  we  consider  the  case  in  which  arc  costs  and  capacities 
are  dependent.  The  random  vectors  Z3  =  (Xj,  Yj)  take  on  values  (cj(  1),  u,(l)),  . . ., 
(c_,(nJ),uJ(nJ))  with  respective  probabilities  p,(l),. . .  ,pj(n,j).  We  further  assume 
that  Uj(l)  >  ...  >  iij(rij)  for  all  j.  Let  Z  =  (Z\, . . . ,  Za).  For  a  fixed  realization  2 
of  Z,  let  C(z)  be  the  cost  of  a  MCF. 

Problem  MC2  Let  d  >  0.  Assuming  the  components  of  Z  are  independent,  eval¬ 
uate  the  probability  P{C( Z)  <  d}  that  there  is  a  flow  of  /  units  whose  cost 


does  not  exceed  d. 


Here  again,  feasible  states  are  those  having  MCF  cost  no  greater  than  d.  Because 
of  the  reverse  ordering  of  the  capacity  variables,  we  have  low  costs  paired  with  high 
capacities  so  that  Zj  is  lower  feasible  for  all  j.  We  will  address  this  problem  in 
Section  4.3 

The  following  result  shows  that  these  problems  are  #P-hard.  It  is  stated  for 
problem  MCI,  but  clearly  implies  that  problem  MC2  is  also  #P-hard. 

Proposition  4.1  The  evaluation  o/P{C'(X|Y)  <  d)  is  #P-hard. 

Proof:  The  evaluation  of  /^{maximum  s  —  t  flow>  d}  is  #P-hard  for  arbitrary  d  >  0 
[3].  Consider  an  instance  of  this  problem.  Put  in  a  return  arc  (t,s)  with  infinite 
capacity.  Let  the  costs  of  all  original  arcs  be  0  and  the  cost  of  arc  ( t,s )  be  —  1. 
Clearly,  a  minimum  cost  flow  in  the  resulting  network  will  ship  as  much  as  possible 
on  arc  ( t,s ),  and  will  thus  yield  a  maximum  s  —  t  flow  in  the  original  network.  In 
fact,  the  events  {C(X|Y)  <  —  d)  and  {maximum  s  —  t  flow>  d}  are  equal. 

Since  the  evaluation  of  ^{maximum  s  —  t  flow  >  d]  is  #P-hard,  so  must  be  the 
evaluation  of  P{C(X|Y)  <  d}.  □ 

In  Section  4.2  we  shall  present  eight  different  GSD  algorithms  for  solving  problem 
MCI,  one  “hybrid”  partitioning  routine,  and  a  GSD  modification  which  simultane¬ 
ously  computes  the  entire  distribution  C(X|Y).  Section  4.3  contains  analogous 
routines  applied  to  problem  MC2  and  the  evaluation  of  P{C( Z)  <  d}.  Numerical 
examples  are  found  in  Section  4.4  and  Section  4.5  offers  concluding  remarks. 

4.2  The  Case  of  Independent  Costs  and  Capac¬ 
ities 

We  consider  here  the  flow  network  whose  arc  costs  X3  and  arc  capacities  Y  j  are 
mutually  independent  random  variables.  The  system  operates  feasibly  if  demand  is 
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satisfied  and  the  constraint  C(X|Y)  <  d  is  satisfied.  We  wish  to  compute  PF  — 

P{C(X|Y)  <  d}. 

In  the  following  subsections,  we  propose  eleven  algorithms  for  solving  problem 
MCI.  We  begin  in  Sections  4.2.1  and  4.2.2  by  describing  two  different  approaches 
to  deriving  feasible  and  infeasible  states.  These  operations,  of  course,  are  the  heart 
of  any  partitioning  routine.  Based  on  each  of  these  approaches,  we  present  in  Sec¬ 
tion  4.2.3  two  GSD  feasible  decompositions  and  two  GSD  infeasible  decompositions, 
and  in  Section  4.2.4  one  hybrid  routine  (in  which  each  interval  is  partitioned  into 
one  arbitrary  feasible  interval,  one  arbitrary  infeasible  interval  and  up  to  4a  —  1 
undetermined  intervals).  In  Section  4.2.5  we  describe  the  simultaneous  evaluation 
of  the  entire  distribution  of  C(X|Y). 

4.2.1  Decomposition  Algorithms  1 

In  this  section  we  describe  techniques  for  deriving  feasible  and  infeasible  cutoff 
states  for  use  in  partitioning  algorithms  for  evaluating  PF .  These  algorithms,  which 
we  collectively  call  the  Minimum  Cost  Flow  Decomposition  Algorithms  I  (MCDA- 
I),  have  the  common  feature  that  feasible  sets  are  removed  by  finding  a  feasible 
state  (k\l)F  derived  from  the  extreme  feasible  point  (a|/?)  and  infeasible  sets  are 
removed  by  deriving  an  infeasible  state  (k]!)1  from  (/?|c*).  We  describe  here  the 
mechanics  for  deriving  ( k\l)F  and  (k\l)1  for  an  undetermined  interval  U  =  [a,/?] 
for  which  C(a\0)  <  d  and  C(0\a)  >  d.  Note  that  a  and  0  are  of  the  form  a  = 
•••,«”)  and  0  =  (/?f, . . . ,  0‘\0^ 

MCDA-I  Feasible  State  Derivation 

Here  we  seek  to  derive  a  feasible  state  (ib|/)F,  starting  with  (a|/?),  that  will  cut 
off  as  large  a  feasible  set  as  is  practical  using  Proposition  2.1. 

We  begin  by  finding  a  minimum  cost  feasible  flow  f°  at  {o\0).  Clearly  if  f°(j)  =  0 
for  some  j  E  A,  we  may  raise  the  cost  of  arc  j  to  Cj(0})  and  reduce  the  capacity  of 
arc  j  to  u}{otj)  without  impacting  optimality.  Also,  if  0  <  f°(j)  <  u3(0j)  we  may 
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reduce  the  capacity  of  arc  j  as  long  as  it  does  not  fall  below 

Having  made  these  simple  changes,  we  are  left  to  decide  whether  or  not  to 
attempt  further  modifications  to  ( fc|/)F.  Of  course,  we  are  free  to  use  (k\l)p  as  it  now 
stands  since  by  Proposition  2.1  any  feasible  point  will  yield  a  feasible  set.  However, 
it  is  likely  that  most  arcs  carry  some  flow,  so  that  in  some  cases  our  previous  changes 
to  (fcj/)F  accomplish  little.  Thus  we  are  tempted  to  try  further  cost  increases  and/or 
capacity  reductions.  At  issue  is  the  fact  that  such  attempts  must  now  come  at  the 
expense  of  re-optimizing.  However,  the  existence  of  efficient  codes  for  performing 
sensitivity  analysis  on  network  flow  problems  eases  our  concerns. 

Since  capacities  have  already  been  addressed  to  a  degree,  we  concentrate  here 
on  raising  arc  costs.  Specifically,  we  attempt  to  raise  as  many  arc  costs  as  possible 
while  remaining  feasible.  Intuitively,  increasing  the  cost  of  a  little-used  arc  should 
have  a  smaller  impact  on  overall  cost  than  increasing  the  cost  of  a  saturated  arc 
(thus  likely  allowing  other  arc  costs  to  be  raised).  For  this  reason,  we  file  the  arcs 
for  which  ac-  <  0C-  in  a  list  ranked  in  increasing  order  of  utilization  (given  initially 
“  f°{j)/u 

In  each  step  we  remove  the  top  arc  from  the  list,  increase  its  cost  one  level, 
reoptimize,  and  verify  that  the  resulting  optimal  flow  still  has  feasible  cost.  If  not, 
the  cost  is  returned  to  the  previous  level,  k,  and  this  arc  is  no  longer  considered.  If 
the  cost  is  feasible,  the  increase  is  retained.  At  this  point,  we  check  if  the  current 
arc  has  zero  flow.  If  so,  its  cost  is  raised  and  its  capacity  is  lowered  as  much  as 
possible.  Otherwise,  we  lower  the  capacity  of  the  arc  being  considered  until  further 
reduction  would  impact  that  arc’s  current  optimal  flow.  If  after  these  changes  the 
arc’s  cost  can  still  be  raised,  the  arc  is  filed  back  in  the  list  (again  in  increasing  order 
of  utilization  of  capacity,  now  at  the  current  (k\l)F)  for  a  subsequent  pass.  We  then 
move  on  to  the  next  arc  to  be  considered.  This  process  continues  until  the  list  is 
empty  or  until  we  achieve  a  cost  which  is  within  a  specified  tolerance  of  the  target 
cost  d. 
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At  that  point  we  have  the  feasible  point  ( k\l)F ,  which  defines  the  feasible  set 
F  =  {(Ar|/)  €  U  :  oj  <  kj  <  kj  for  j  E  A 

h  <  h  <  fij  for  j  6  A}.  (4.1) 

This  feasible  state  derivation  will  be  used  in  three  different  partitioning  routines, 
as  we  shall  describe  in  Sections  4.2.3  and  4.2.4. 

MCDA-I  Infeasible  State  Derivation 

Here  we  seek  to  derive  an  infeasible  state  (fc|/);,  starting  with  the  infeasible  state 
(/?|a),  that  will  cut  off  as  large  an  infeasible  set  as  is  practical  using  Proposition  2.1. 

We  begin  by  finding  a  minimum  cost  flow  f°  (with  infeasible  cost)  at  (/?| o). 
All  arcs  j  for  *vuich  aj  <  /?|  or  a“  <  (i.e.,  all  arcs  whose  cost  or  capacity  can 

be  changed)  are  filed  in  a  list  and  are  considered  in  increasing  order  of  utilization, 
given  as  /°(j)/uJ(a“).  The  reasoning  behind  this  is  that  in  this  way  we  we  first 
look  at  changing  the  parameters  ofthe  least-used  arcs.  (Recall  that  we  would  like 
for  the  infeasible  point  to  agree  with  (a\0)  in  as  many  positions  as  possible  so  as  to 
minimize  the  number  of  the  resulting  undetermined  sets.) 

At  each  step  we  remove  the  top  arc  in  the  list  and  attempt  to  lower  its  cost  and 
raise  its  capacity  one  level  each.  After  each  of  these  changes,  we  reoptimize  and 
check  if  the  resulting  state  is  still  infeasible.  If  not,  the  previous  cost  or  capacity  is 
restored.  But  if  either  of  these  changes  can  be  retained,  and  either  the  cost  or  the 
capacity  can  still  be  changed  (i.e.,  cost  has  not  yet  been  lowered  to  the  level  ac  or 
the  capacity  has  not  been  raised  to  the  level  /?"),  the  arc  will  be  filed  back  in  the  list 
for  future  consideration,  in  increasing  order  of  its  current  utilization.  This  process 
continues  until  no  arcs  remain  to  be  considered  or  the  MCF  cost  drops  to  within  a 
specified  tolerance  of  d. 

At  that  point,  we  have  the  infeasible  state  ( fcj/)7 ,  which  defines  the  infeasible  set 

/  =  {(Ar|/)  €  U  :  k}  <  kj  <  fl'  for  j  €  A 

o'-  <lj  <  lj  for  >€.4}.  (4.2) 
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This  infeasible  state  derivation  will  be  used  in  three  different  partitioning  rou¬ 
tines,  as  we  shall  see  after  presenting  alternative  feasible  and  infeasible  state  deriva¬ 
tions  in  Section  4.2.2. 

4.2.2  Decomposition  Algorithms  II 

In  this  section  we  discuss  an  alternative  approach  to  deriving  feasible  and  infea¬ 
sible  cutoff  states  for  use  in  partitioning  algorithms  for  evaluating  PF ,  referred  to 
collectively  as  the  Minimum  Cost  Flow  Decomposition  Algorithms  II  (MCDA-II). 
The  feasible  and  infeasible  state  derivations  given  below  have  the  feature  that  fea¬ 
sible  states  are  derived  starting  at  the  infeasible  extreme  (0\ a)  and  infeasible  states 
are  derived  from  the  feasible  extreme  (»|/3).  We  describe  here  the  procedures  for 
deriving  feasible  and  infeasible  state  points  in  an  undetermined  interval  U  —  [a,  f3] 
for  which  C(a|/3)  <  d  and  C(fl\a)  >  d.  Note  that,  again,  a  and  /3  are  of  the  form 
used  in  Section  4.2.1. 

MCDA-II  Feasible  State  Derivation 

Here  we  seek  to  determine  a  feasible  point  ( k\l)F  by  moving  from  the  infeasible 
extreme  (/?| a). 

We  begin  by  finding  /°,  a  minimum  cost  flow  at  (/ 3\a ).  Note  that  for  reasons 
presented  in  Section  2.5,  we  wish  to  move  away  from  (f3 |a)  in  as  few  arcs  as  possible. 
Therefore  we  shall  consider  arcs  in  decreasing  order  of  utilization  (defined  initially 
35  /°(j)/uj(at“)).  The  reasoning  here  is  that  if  we  first  raise  capacities  and  lower 
costs  for  heavily-used  arcs,  we  might  achieve  a  feasible  cost  quickly  and  leave  many 
arcs  set  at  the  (/ 3\a )  levels.  Of  course,  as  soon  as  a  feasible  cost  is  achieved  we  stop 
since  our  goal  is  then  accomplished. 

Thus  we  first  identify  the  most  heavily-used  arc.  If  the  cost  of  this  arc  can  be 
reduced,  we  reduce  it  all  the  way,  reoptimize  and  compute  a  new  MCF  cost.  If 
the  cost  is  feasible,  we  stop.  Otherwise,  we  determine  whether  or  not  this  arc  is 
still  most  heavily-used.  If  it  is,  we  attempt  to  raise  is  capacity  all  the  way,  again 


81 


checking  for  the  effect  on  the  MCF  cost,  If  the  current  arc  is  no  longer  the  most 
attractive  arc,  we  begin  anew  lowering  the  cost  of  the  most  attractive  arc,  and  the 
cycle  is  repeated. 

Two  practical  points  bear  mentioning  here.  The  first  point  is  that  only  arcs 
with  positive  flow  are  considered  for  the  changes  described  above.  As  a  result,  the 
method  as  stated  can  stall.  That  is,  we  can  arrive  at  a  state  which  still  has  infeasible 
cost,  and  for  which  all  utilized  arcs  have  already  been  moved  as  much  as  possible 
from  (/?|a).  In  these  situations  we  lower  the  costs  of  unused  arcs  as  much  as  possible 
until  the  situation  is  remedied  (which  it  must  eventually  be  since  we  know  (a|/?)  is 
feasible).  We  first  consider  such  arcs  with  lowest  possible  costs  and  continue  in  this 
order  until  the  main  method  can  continue. 

The  second  point  is  that  the  determination  of  most-used  arcs  and/or  lowest- 
cost  unused  arcs  is  done  with  almost  no  additional  effort.  This  is  so  because  we 
determine  this  information  in  the  course  of  looping  through  the  arcs  to  compute  the 
new  MCF  cost  after  reoptimizing.  Thus,  these  arcs  are  always  at  hand  whenever 
they  are  needed. 

This  feasible  state  determination  will  be  used  in  three  different  partitioning  rou¬ 
tines,  as  we  shall  see  in  Sections  4.2.3  and  4.2.4. 

MCDA-II  Infeasible  State  Derivation 

Here  we  seek  to  derive  an  infeasible  point  {k\I)1  by  starting  at  the  feasible  extreme 
(<*(/?)  and  changing  costs  and/or  capacities  of  only  as  many  arcs  as  is  necessary  to 
achieve  a  MCF  cost  exceeding  d.  We  wish  for  (fc|/);  to  agree  with  (a|/3)  in  as  many 
positions  as  possible  so  as  to  minimize  the  number  of  the  resulting  undetermined 
sets. 

We  begin  by  finding  /°,  a  minimum  cost  feasible  flow  at  (a|/3).  As  in  MCDA-II 
feasible  decomposition,  we  identify  the  most  heavily-used  arc  at  the  current  best 
point  by  raising  costs  and/or  lowering  capacities  in  an  effort  to  achieve  an  infeasible 
cost.  We  stop  as  soon  as  a  cost  greater  than  d  is  achieved.  The  same  practical 
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concerns  mentioned  above  (stalling,  work  done  to  find  new  arcs)  are  also  relevant 
here. 

This  infeasible  state  derivation  will  be  used  in  three  different  partitioning  rou¬ 
tines,  as  we  now  discuss. 

4.2.3  GSD  Routines  for  the  Independent  Case 

In  Sections  4.2.1  and  4.2.2  we  presented  methods  for  deriving  feasible  and  infea¬ 
sible  cutoff  states.  These  point  definitions  were  made  for  the  purpose  of  removing 
feasible  and  infeasible  intervals  in  accordance  with  Proposition  2.1.  We  describe 
here  GSD  routines  that  will  use  these  feasible  and  infeasible  point  definitions  to 
evaluate  PF  =  P{C(X|Y)  <  d}.  The  efficiency  of  these  routines  will  be  tested 
Section  4.4.  In  each  case  below  we  describe  the  partitioning  of  an  undetermined 
interval  U  =  [a,  /?]  for  which  C(a\f3)  <  d  and  C(/?|a)  >  d. 

GSD  Feasible  Decompositions 

The  most  straightforward  feasible  decomposition  consists  of  simply  deriving  a 
feasible  point  {k\l)F,  removing  the  resulting  feasible  set  defined  by  (4.1),  and  par¬ 
titioning  the  remainder  of  U  into  at  most  2a  new  undetermined  intervals.  For  each 
arc  j  we  have  two  undetermined  intervals, 

Uj  =  { ( A:  |  / )  E  U  :  ft  <  k,  <  kt  for  i  <  j 

kj  +  l  <kj<  /i;  (4.3) 

atc  <  k,  <  (if  for  i  >  j 

«“</,<  ft  for  1  €  A} 

and 

Uj  —  {(fc|f)  G  U  :  ctl  <  ki  <  ki  for  i  E  A 

h  <  l,  <  ft  for  i  <  j 

a?<h<ij-  l  (4.4) 
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< *?<h<0 t  for  i  >  jj. 


But  in  many  cases  we  wish  to  also  remove  infeasible  sets  in  each  iteration.  Section  2.5 
suggests  that  this  is  best  accomplished  by  finding  limiting  feasible  indices  k3  and  l3 
as  given  in  Definition  2.8.  Note  that  some  of  these  are  already  known.  If  k3  =  fic- 


for  any  j  then  kj  =  0'.  Similarly,  —  o'-  if  /,  =  aj. 

After  we  find  all  2 a  indices, 

the  standard  GSD  feasible  decomposition  of  Section  2.4  removes  up  to  2 a  infeasible 

sets  of  the  form 

/;  =  {(*|/)  e  u : 

<  h  <  ki 

for  i  <  j 

h  +  1  <  kj  <  /3cj 

(4.5) 

<  <  h  <  0\ 

for  i  >  j 

<*“</,<  0? 

for  i  €  A} 

and 

I?  =  {(k\l)eU: 

act  <  kj  <  k, 

for  i  e  A 

U  <  U  <  # 

for  i  <  j 

<b<h-l 

(4.6) 

af  <  U  <  0? 

for  i  >  j  } 

and  up  to  2a  new  undetermined  intervals  given  by 

U;  =  {(*10  €  U  : 

5;  k,  <  fc, 

for  i  <  j 

kj  -(-  1  ^  k3  ^  kj 

(4.7' 

ai  <  ki  <  k{ 

for  i  >  j 

/,</,<  0t 

for  i  £  A} 

and 

u;  =  im  e  u : 

act  <  ki  <  k{ 

for  i  £  A 

ii<h<  0t 

for  i  <  j 
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h<b<h- 1 

h  <U<  ft  for  i  >  j } . 


(4.8) 


Thus,  using  MCDA-I  or  MCDA-II  feasible  point  derivations  with  or  without  limiting 
feasible  indices  gives  us  four  GSD  feasible  routines,  each  of  which  follows  the  basic 
form  of  what  was  called  Pure  Option  1  in  Section  2.4. 


GSD  Infeasible  Decompositions 


As  with  feasible  decomposition,  we  have  the  option  of  using  limiting  indices  or 
not.  If  we  decline  to  do  so,  then  we  simply  derive  the  infeasible  state  ( k |/);  by 
one  of  the  two  methods  given  above,  remove  the  infeasible  set  /  given  in  (4.2),  and 
partition  the  remainder  of  U  into  up  to  2a  undetermined  intervals  g.  'en  by 


and 


U]  =  {(*10  e  U  :  hi  <  ki  <  ft  for  i  <  j 

aCj  <  kj  <  kj  —  1 

<  ki  <  ft  for  i  >  j 

«“</,<  /?“  for  i  G  A } 


(4.9) 


Uf  =  {(*|/)  €  U  :  h<ki<P‘  ioviZA 

a'-  <li  <  li  for  i  <  j 

h + 1  <  h  <  pi 

<*i<li<Pi  for  *•>;}. 


(4.10) 


If  we  choose  to  also  remove  feasible  sets  during  each  iteration,  we  first  determine 
limiting  infeasible  indices  kj  and  lj  as  in  Definition  2.9.  (Again,  if  kj  =  ac}  we  may 
set  kj  =  aCj  and  if  J}  =  PJ  we  may  set  lj  =  PJ.)  Once  all  indices  are  determined,  we 
may  remove  up  to  2a  feasible  sets  of  the  form 

Fcj  =  {(Jfc|0  €  U  :  ki  <  ki  <  P\  for  i  <  j 


85 
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and 


aCj  <kj  <kj  —  1 

a-  <  ki  <  P-  for  i  >  j 

a*  <U<  ft  for  i  €  >1} 

FJ  =  {(k\l)  eU:  k<k<ft  forieA 
atu  <h<  /j  for  t  <  j 

h  + 1  <  h  <  ft 
ft  <U<ft  ^  *  >  >} 


(4.11) 


(4.12) 


and  partition  the  remainder  of  U  into  up  to  2a  new  undetermined  sets  given  by 


U]  =  {(*|/)  €  U  :  *,  <  ki  <  ft 

for  i  <  j 

kj  <  kj  <kj  —  1 

(4.13) 

ki  <  ki  <  ft 

for  t  >  j 

ft  <  u  <  h 

for  i  (E  .4} 

U»  =  {(*|0  <=  U  :  h<ki<Pt 

for  i  e  A 

<u<  u 

for  i  <  j 

h  + 1  <  ij  <  h 

(4.14) 

ft  <  /,  <  k 

for  t  >  j}. 

Thus  we  have  four  different  GSD  infeasible  decomposition  routines,  each  follow¬ 
ing  the  general  description  of  Pure  Option  2  in  Section  2.4. 

4.2.4  Hybrid  Decomposition  Routine  for  the  Independent 
Case 


In  Section  2.5  we  were  cautioned  that  deviating  from  GSD  would  most  likely 
restilt  in  diminished  performance  since  such  a  departure  was  proved  to  cause  an 


» 


I 


I 


ft 


I 


I 


» 


» 


86 


I 


ft 


unnecessary  proliferation  of  undetermined  intervals.  Nevertheless,  since  many  other 
factors  play  a  part  in  determining  a  routine’s  ultimate  effectiveness,  it  is  worth¬ 
while  to  attempt  other  kinds  of  partitioning  schemes.  In  that  spirit,  we  describe 
here  a  routine  that  is  intuitively  very  appealing,  but  that  will  file  up  to  4a  —  1  new 
undetermined  intervals  in  each  iteration  (instead  of  up  to  2 a  as  in  GSD  routines). 
This  procedure,  which  we  have  called  a  hybrid  routine,  removes  the  largest  prac¬ 
tical  feasible  interval  and  the  largest  practical  infeasible  interval  in  each  iteration. 
Specifically,  in  this  routine  we  determine  a  feasible  state  (k\l)F  by  either  MCDA-I  or 
MCDA-II  feasible  state  definition  and  an  infeasible  state  (k\l)!  by  the  corresponding 
infeasible  state  definition.  We  remove  the  resulting  feasible  set  F  as  given  in  (4.1), 
adding  its  probability  to  the  accumulated  lower  bound  for  PF ,  and  the  infeasible 
set  I  as  given  in  (4.2),  and  decompose  the  remainder  of  U  in  the  manner  described 
in  the  proof  of  Proposition  2.3. 

Thus  we  have  two  different  hybrid  routines  (one  corresponding  to  MCDA-1  and 
one  corresponding  to  MCDA-II). 

We  now  conclude  our  discussion  of  the  case  of  independent  arcs  by  describing 
the  algorithm  MCDIST. 

4.2.5  Computing  the  Entire  Distribution  of  MCF  Costs 

The  algorithms  presented  thus  far  are  well-suited  to  situations  in  which  a  limit  on 
the  cost  of  an  optimal  flow  exists  a  priori  (e.g.,  budgetary  constraints,  availability  of 
materials,  imposed  specifications).  If  that  is  not  the  case,  it  is  useful  to  completely 
characterize  the  probabilistic  behavior  of  C(X|Y),  including  the  computation  of 
the  mean  and  higher  moments  of  C(X|Y).  In  this  section  we  describe  MCDIST,  a 
routine  that  provides  precisely  this  information  by  producing  the  probability  distri¬ 
bution  of  C(X|Y). 

In  a  slight  departure  from  straight  GSD,  MCDIST  is  not  concerned  with  cal¬ 
culating  a  single  probability  measure.  Rather,  it  is  concerned  with  computing  the 
probability  function  q(d)  =  P{C(X|Y)  =  d),  d>  0.  Lower  bounds  for  these  values, 
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denoted  qi(d),  are  initialized  to  zero  and  accumulate  the  appropriate  probabilities 
as  decomposition  proceeds.  In  the  following  description,  U  consists  of  those  unde¬ 
termined  sets  not  yet  considered.  As  with  the  procedure  STDIST,  we  assume  that 
MCF  costs  will  be  integral  (which  we  enforce  by  requiring  integral  cost  and  capacity 
values). 

Procedure  MCDIST 

Step  1  Find  C/„  =  C(l, . . . ,  l|ni, . . .  ,na)  and  Chi  =  C(mi, . . .  ,m0|l, . . . ,  1).  For 
d  =  Cio, . . . ,  Chi ,  set  qi(d)  =  0.  Set  U  =  0  and  U  =  Cl. 

Step  2  Let  a  and  0  denote,  respectively,  the  lower  and  upper  endpoints  of  U.  Find 
/“,  a  minimum  cost  flow  at  (a|/?),  with  cost  C{a\0). 

Step  3  Let  (k\l)  be  derived  from  (a\0)  by  pushing  costs  and  capacities  of  unused 
arcs  to  the  opposite  extreme  and  reducing  unused  capacities. 

Step  4  Set  qi(C(a\0))  =  qt(C{oc\0))  +  P{F},  where  F  is  as  defined  by  (4.1). 

Step  5  For  j  €  A:  If  kj  /  0^,  file  the  interval  Uj  defined  by  (4.3)  in  U.  If  ^  a“, 
file  the  interval  UJ  defined  by  (4.4)  in  U. 

Step  6  If  U  =  0,  stop.  Otherwise,  remove  a  set  U  from  U  and  go  to  step  2. 

Note  that  in  step  3  there  really  isn’t  room  for  improving  (k\J)  (i.e.,  for  pushing  it 
further  from  (a|/?))  since  increasing  the  cost  of  an  arc  which  carries  flow  necessarily 
changes  the  cost  of  a  minimum  cost  flow. 

As  usual,  MCDIST  may  be  terminated  prior  to  complete  partitioning  of  fi  to 
produce  bounds.  One  can  stop,  for  example,  when  the  bounds  are  sufficiently  close 
or  when  a  predetermined  number  of  sets  have  been  decomposed.  To  produce  bounds 
on  F(d)  =  P{C(X|Y)  <  d}  for  all  d,  we  substitute  step  6  for  step  6,  and  add  steps 
7  and  8. 


Step  6  If  U  —  0,  construct  the  cumulative  distribution  function,  and  stop.  If  a 
stopping  condition  is  met,  go  to  step  7.  Otherwise,  remove  a  set  V  from  14 
and  go  to  step  2. 

Step  7  For  d  Clo, . . . ,  C\t:  let  Ft(d)  =  FM(d)  =  ?/(*)• 

Step  8  For  all  £/  €  £7:  Let  cv  and  0  denote,  respectively,  the  lower  and  upper 
endpoints  of  U.  Find  C( a\0)  and  C(0\a). 

For  C{a\0)  <d<  C{0\a):  Set  Fu(d)  =  Fu(d)+  P{U}. 

For  C{0\a)  <d<  Chi:  Set  Fu(rf)  =  Fu(d)  +  P{U)  and  F,{d)  =  Ft{d)  +  P{U}. 

The  performance  of  this  algorithm  will  be  evaluated  in  Section  4.4. 

4.3  The  Case  of  Dependent  Costs  and  Capacities 

Thus  far  we  have  considered  minimum  cost  flow  problems  in  which  arc  costs  ai.d 
capacities  are  given  as  mutually  independent  random  variables.  In  practice,  though, 
it  is  sometimes  the  case  that  the  cost  and  upper  bound  for  a  given  arc  vary  together. 

Example  4.1  A  link  in  a  transportation  network  can  move  20  units  of  commodity 
per  unit  time  at  an  average  cost  of  $.15  per  unit  when  it  is  in  perfect  condition. 
The  link  can  sustain  minor  damage  and  still  operate  with  diminished  capacity  of 
15  units  and  a  moderate  increase  in  cost  to  $.20  (where  the  increased  cost  is  due 
to  additional  required  monitoring  of  the  link).  Serious  damage  to  the  link  causes 
its  capacity  to  drop  to  only  5  units  and  increases  the  average  cost  per  unit  to  $.30 
since  at  that  point,  a  worker  is  dispatched  to  make  repairs  on  the  link.  Finally,  the 
link  can  be  completely  failed  (its  capacity  is  reduced  to  0  units). 

In  order  to  model  situations  such  as  that  given  in  Example  4.1,  we  associate  with 
each  arc  a  vector-valued  random  variable  Zj  =  (Xj,  Yj),  where  Xj  and  Y}  represent 
cost  and  upper  bound,  repectively,  as  in  Section  4.2,  but  where  Xj  is  a  deterministic 
function  of  Yj.  In  Example  4.1,  Zj  takes  on  values  (.15,20),  (.20,15),  (.30,5)  and 
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(*,0).  (The  cost  per  unit  associated  with  capacity  0  is  immaterial.)  As  mentioned 
earlier,  we  require  that  c ^ ( 1 )  <  •••  <  e;(n;)  and  »!>(1)  >  •••  >  Uj(jiy)  so  that  we 
may  classify  Zj  as  a  lower  feasible  variable.  The  state  space  ft  of  Z  =  (Z j,. . . ,  Za) 
consists  of  ntsi  nj  state  points  of  the  form  z  —  ((c,(/ t),«i(/i)), ••  • , (c0(/a), ua{la))), 
where  l,  6  {1,...  ,»»,■}  for  each  j  G  A.  As  before,  we  shall  refer  to  this  state  point 
as  /  —  (/j, . 

In  the  remainder  of  this  section,  we  present  a  number  of  algorithms  for  solving 
problem  MC2.  To  the  extent  possible,  they  will  largely  be  adaptations  of  the  routines 
given  in  Section  4.2  for  the  case  of  independent  costs  and  capacities.  However,  we 
shall  see  that  some  flexibility  is  lost  in  translating  those  ideas  into  the  present 
setting.  We  will  now  be  unable  to  change  only  the  cost  or  only  the  capacity  of 
an  arc.  For  example,  we  cannot  reduce  unused  arc  capacity  without  impacting  an 
optimal  solution  since  that  arc’s  cost  will  necessarily  increase. 

In  Section  4.3.1  we  give  dependent  versions  of  feasible  and  infeasible  state  deriva¬ 
tions  corresponding  to  MCDA-I.  Section  4.3.2  gives  derivations  for  the  dependent 
MCDA-Il.  GSD  routines  for  solving  the  present  problem  are  given  in  Section  4.3.3, 
and  “hybrid”  algorithms  are  given  in  Section  4.3.4.  In  Section  4.3.5  we  present  the 
dependent  version  of  MCDIST  for  computing  the  entire  distribution  of  C( Z).  The 
performance  of  all  these  routines  will  be  evaluated  in  Section  4.4. 

4.3.1  Decomposition  Algorithms  I 

We  present  here  feasible  and  infeasible  state  derivations  to  be  used  in  a  number 
of  partitioning  algorithms  for  solving  problem  MC2.  W’e  seek  to  partition  an  unde¬ 
termined  interval  U  —  [a,/?]  for  which  C(a)  <  d  and  C{0)  >  d.  These  derivations 
have  the  common  feature  that  feasible  cutoff  states  will  be  found  by  starting  at  the 
feasible  extreme  a  (which  corresponds  to  lowest  cost  paired  with  highest  capacity 
for  each  arc)  and  infeasible  cutoff  states  will  be  derived  from  the  infeasible  extreme 
0  (highest  cost  paired  with  lowest  capacity). 
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MCDA-I  Feasible  State  Derivation 

Here  we  wish  to  derive  a  feasible  state  F  by  moving  from  the  extreme  feasible 
point  a.  We  begin  by  finding  /“,  a  minimum  cost  feasible  flow  at  a.  As  with  the 
independent  MCDA-I,  we  may  simultaneously  raise  the  arc  cost  and  lower  the  arc 
capacity  of  any  arc  that  does  not  carry  flow  in  f°  without  impacting  optimality. 
(Recall  that  cost  and  capacity  form  a  dependent  pair  so  that  this  operation  consists 
of  moving  to  higher- indexed  state  vectors.)  After  doing  this,  we  file  the  remaining 
arcs  j  (for  which  atj  <  fy)  in  a  heap  sorted  in  increasing  order  of  utilization  (given 
as  f°(j)/uj(aj)).  As  explained  before,  we  do  this  so  as  to  make  If  close  to  0j  for  as 
many  arcs  as  possible  to  cut  down  on  the  number  of  undetermined  sets  filed. 

We  consider  arcs  in  this  order  and  attempt  to  simultaneously  increase  arc  cost 
and  reduce  arc  capacity  a  single  level.  After  each  such  move,  we  re-optimize  and 
check  the  new  MCF  cost.  If  the  move  has  caused  the  cost  to  exceed  d,  we  restore  the 
previous  cost  and  capacity  level  and  eliminate  the  arc  from  further  consideration. 
Otherwise,  the  change  is  retained  and  we  re-file  the  arc  (as  long  as  we  have  not 
raised  its  cost  and  lowered  its  capacity  as  much  as  possible)  for  a  future  pass.  This 
process  continues  until  no  arcs  remain  to  be  considered  or  the  MCF  cost  is  raised 
to  within  a  specified  tolerance  of  the  target  cost  d. 

MCDA-I  Infeasible  State  Derivation 

Here  we  seek  to  derive  an  infeasible  point  V  by  moving  from  the  infeasible 
extreme  /?.  Our  goal  is  to  use  this  point  to  cut  off  as  large  a  feasible  set  as  possible 
without  causing  unnecessary  proliferation  of  new  undetermined  sets.  As  such,  we 
will  again  try  to  move  as  many  arcs  as  possible  toward  the  level  o  by  first  considering 
little-used  arcs  and  proceeding  in  order  of  increased  utilization. 

We  begin  by  finding  /°,  a  minimum  cost  flow  (with  cost  exceeding  d)  at  the  point 
/?.  We  then  consider  arcs  in  increasing  order  of  utilization.  For  each  arc  considered, 
we  lower  its  cost  and  raise  its  capacity  a  single  level,  keeping  this  change  only  if  the 
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resulting  MCF  still  has  infeasible  cost.  If  that  is  the  case,  and  there  is  still  room 
to  move  this  arc’s  cost  and  capacity,  we  file  it  back  in  a  list  for  future  passes.  This 
process  continues  until  there  are  no  more  arcs  to  consider  or  the  MCF  cost  falls  to 
within  a  specified  tolerance  of  d. 

4.3.2  Decomposition  Algorithms  II 

This  section  contains  feasible  and  infeasible  state  derivations  to  be  used  in  par¬ 
titioning  algorithms  for  solving  problem  MC2.  The  descriptions  below  give  details 
for  finding  feasible  and  infeasible  cutoff  states  for  an  undetermined  set  U  =  [a,  0} 
having  (7(a)  <  d  and  C{0)  >  d.  Feasible  states  will  be  derived  starting  at  the 
infeasible  extreme  0  and  infeasible  states  will  be  derived  from  a.  The  motivation 
for  this  approach  is  to  produce  feasible  cutoff  states  that  agree  with  0  in  as  many 
positions  as  possible  and  infeasible  cutoff  states  close  to  a.  As  such,  since  we  wish  to 
move  few  arcs  away  from  their  starting  extreme,  we  shall  consider  arcs  in  decreasing 
order  of  utilization. 

MCDA-II  Feasible  State  Derivation 

Here  we  derive  a  feasible  point  lF  from  the  infeasible  extreme  0.  We  begin  by 
finding  /°,  a  minimum  cost  flow  (with  cost  exceeding  d)  at  the  point  0.  We  wish 
to  look  first  at  heavily-used  arcs.  As  in  the  independent  version  of  MCDA-II,  at 
all  times  we  have  available  the  current  most  utilized  arc.  We  reduce  the  cost  and 
increase  the  capacity  of  this  arc  all  the  way,  and  then  re-optimize  to  determine  the 
effect  of  this  change  on  the  MCF  cost.  We  continue  in  this  fashion  until  reoptimiza¬ 
tion  reveals  that  a  feasible  cost  is  achieved.  At  that  point  we  stop. 

MCDA-II  Infeasible  State  Derivation 

We  seek  to  derive  an  infeasible  state  point  l1  by  starting  at  the  feasible  extreme 
a.  The  method  by  which  we  do  this  is  essentially  the  same  as  that  given  above  for 
feasible  states.  Specifically,  we  find  a  minimum  cost  feasible  flow  f°  at  a  and  then 
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move  away  from  this  point  in  as  few  arcs  as  possible.  We  accomplish  this  by  always 
considering  the  most  heavily-used  arc.  We  stop  as  soon  as  a  MCF  cost  exceeding  d 
has  been  achieved. 

4.3.3  GSD  Routines  for  the  Dependent  Case 

In  Sections  4.3.1  and  4.3.2  we  described  techniques  for  deriving  feasible  and 
infeasible  states  for  the  minimum  cost  flow  problem  with  dependent  arc  costs  and 
capacities.  These  point  definitions  were  made  for  the  purpose  of  removing  feasible 
and  infeasible  intervals  in  partitioning  algorithms  for  solving  problem  MC2.  We 
describe  here  GSD  routines  for  this  problem.  The  examples  in  Section  4.4  will  test 
the  performance  of  these  routines.  In  each  case  below  we  detail  the  partitioning  of 
an  undetermined  interval  U  =  [«,  0]  for  which  C(a)  <  d  and  C(P)  >  d. 

GSD  Feasible  Decompositions 

We  describe  first  the  simplest  GSD  feasible  decomposition,  in  which  either 
MCDA-I  or  MCDA-II  is  used  to  derive  a  feasible  state  lF ,  the  feasible  interval 

F  =  {l£U:  forje.4}  (4.15) 

is  removed,  and  the  remainder  of  U  is  split  into  up  to  a  undetermined  intervals  given 

by 

Uj  =  {/  G  U  :  a,  <  U  <  li  for  i  <  j 

h  +  1  <  h  <  Pi  (4.16) 

a,  <  /,  <  Pi  for  i  >  j}. 

However,  since  the  dependence  of  cost  and  capacity  limit  the  depth  to  which 
we  may  push  a  feasible  state,  we  generally  choose  to  also  remove  infeasible  sets  by 
defining  limiting  feasible  indices  \j  as  given  in  Definition  2.8.  With  the  o  indices 
identified,  we  remove  F  as  above  and  up  to  a  infeasible  intervals 

Ij  =  {l  €  U  :  a,  <  li  <  li  for  i  <  j 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


ft 
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h  +  1  <  h  <  Pi  (4.17) 

<Xi<U<Pi  for  i  >  j}, 

partitioning  the  remainder  into  up  to  a  undetermined  sets  given  by 

Uj  =  {l  e  U  :  ati  <  U  <  for  k  <  j 

lj  +  1  <  lj  <  lj  (4.18) 

a,  <  /,  <  li  for  i  >  j}. 

Thus,  using  either  MCDA-I  or  MCDA-II  with  or  without  limiting  feasible  indices 
yields  four  GSD  routines,  each  following  the  basic  form  of  Pure  Option  1  described 
in  Section  2.4. 

GSD  Infeasible  Decompositions 

As  with  feasible  decomposition,  we  have  the  option  of  using  limiting  infeasible 
indices  or  not.  If  we  decline  to  do  so,  then  we  simply  derive  the  infeasible  state  V 
by  one  of  the  two  methods  presented,  remove  the  infeasible  set 

I={leU:  for;  e  A}  (4.19) 

and  divide  the  remainder  of  U  into  up  to  a  undetermined  sets 

Uj  =  {/  6  U  :  li  <  li  <  ^  for  i  <  j 

Qj  <  lj  <ij-  1  (4.20) 

<*i<li<Pi  for  i  >  j). 

If  we  choose  to  also  remove  feasible  sets  during  each  iteration,  we  first  find 
limiting  infeasible  indices  lj  given  in  Definition  2.9.  In  that  case,  /  is  removed  as 
above  along  with  up  to  a  feasible  intervals 

Fj  =  {l  £  U  :  /,  <U<  Pi  for  *  <  j 

Qj  <lj<ij-  1  (4.21) 

<*i  <1,  <  P,  for  i  >  j), 
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and  up  to  a  undetermined  sets  given  by 

U}  =  {/  €  U  :  U<  U  <  Pi  for  i  <  j 

h  <i,<b- 1  (4.22) 

a,  <  /,  <  k  for  t  >  j}. 

Thus  we  have  four  GSD  infeasible  decomposition  schemes,  each  following  the 
general  description  of  Pure  Option  2  in  Section  2.4. 

4.3.4  Hybrid  Decomposition  Routine  for  the  Dependent 
Case 

As  in  the  case  of  independent  costs  and  capacities,  in  each  iteration  we  attempt  to 
remove  from  the  current  undetermined  interval  U  =  [a,  0}  a  single  arbitrary  feasible 
set  and  a  single  arbitrary  infeasible  set.  The  tradeoff  to  this  intuitive  idea  (removing 
as  much  of  the  feasible  and  infeasible  regions  as  possible)  is  that  the  remainder  of 
U  can  require  as  many  as  2a  —  1  undetermined  intervals  for  its  partition. 

Basically,  this  hybrid  routine  derives  a  feasible  point  lF  by  either  MCDA-I  or 
MCDA-II  feasible  state  derivation  and  an  infeasible  state  ll  by  the  corresponding 
infeasible  derivation.  TLe  probability  of  the  feasible  set  F  (as  given  in  (4.15))  is 
added  to  the  lower  bound  for  PF,  the  infeasible  set  /  (as  in  (4.111)  is  discarded,  and 
U\(F  U  I)  is  partitioned  as  in  the  proof  of  Proposition  2.3. 

4.3.5  Computing  the  Entire  Distribution  of  MCF  Costs 

As  we  indicated  earlier,  it  is  often  useful  to  derive  the  entire  distribution  of  C(Z). 
The  procedure  here  is  completely  analogous  to  that  given  in  Section  4.2.5  with  the 
exception  of  the  derivation  of  our  feasible  point  in  step  3.  The  only  move  we  can 
make  in  pushing  in  from  a  while  keeping  the  MCF  cost  constant  at  C(a)  is  with 
the  unused  arc*,.  We  cannot  decrease  excess  capacity  of  arcs  that  carry  flow  since 
this  would  also  increase  the  arc  cost  and  hence  the  MCF  cost. 
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With  V  thus  defined,  we  remove  the  feasible  set  F  as  in  (4.15)  and  up  to  a 
undetermined  sets  Uj  defined  by  (4.16). 

4.4  Stochastic  MCF  Examples 

In  this  section  we  demonstrate  the  performance  of  the  various  routines  devel¬ 
oped  in  the  present  chapter.  We  test  eight  GSD  routines,  two  hybrid  routines  and 
MCDIST  both  for  the  case  of  independent  arc  costs  and  capacities  and  for  the  case 
of  dependent  pairs  of  arc  costs  and  capacities,  In  each  case,  we  label  GSD  routines 
that  do  not  use  limiting  indices  by  “1”  and  those  using  limiting  indices  by  “2.” 
Our  main  goal  is  to  compare  the  various  approaches  and  evaluate  their  efficiency  in 
producing  bounds.  All  codes  were  implemented  in  FORTRAN  77,  and  were  com¬ 
piled  and  executed  on  a  SUN  SPARCstation  2.  CPU  times  reported  are  in  seconds. 
Minimum  cost  flows  were  found  using  the  relaxation  codes  RELAXT  and  SENSTV, 
developed  by  Bertsekas  and  Tseng  [9].  In  all  applications,  the  list  of  undetermined 
sets  was  kept  as  a  heap  sorted  by  probability  until  the  most  probable  undetermined 
set  had  probability  less  than  10-2°.  In  the  tables  below  we  denote  the  sum  of  the 
probabilities  of  the  undetermined  intervals  remaining  in  the  list  U  (the  difference 
between  upper  and  lower  bounds)  by  P{U). 

Example  4.2  (Example  1.3  revisited)  The  network  in  Table  4.1  has  10  nodes  and 
21  directed  arcs.  Each  arc  has  between  two  and  four  possible  costs  and  between  one 
and  four  possible  capacities.  The  cardinality  of  the  state  space  is  roughly  9  x  1018. 
We  wish  to  ship  15  units  from  node  s  =  1  to  node  t  =  10.  The  range  of  possible 
optimal  flow  costs  is  from  1591  to  3718.  Note  that  this  data  set  differs  from  most 
we  have  used  in  that  it  does  not  have  probability  as  heavily  concentrated  at  the 
status  quo.  Later  in  this  section,  we  re-evaluate  this  problem  with  alternative  input 
distributions.  We  use  the  ten  decomposition  routines  described  above  to  compute 
P{C{X)  <  1900}  and  P{C(X)  <  3500}.  Results  for  F(1900)  are  in  Table  4.3,  for 
F(3500)  in  Table  4.4.  We  see  that  hybrid  routines  did  not  perform  as  well  as  the 


Table  4.1:  Arc  Cost  Distributions  for  Example  4.2. 


Arc 

Cost  (Probability) 

1  =  (1,2) 

70  (0.5) 

73  (0.3) 

94  '0.2) 

2  =  (1,3) 

25  (0.5) 

35  (0.4) 

82  (0.1) 

3  =  (1,4) 

42  (0.5) 

48  (0.3) 

61  (0.2) 

4  =  (2,5) 

26  (0.4) 

31  (0.3) 

55  (0.2) 

88  (0.1) 

5  =  (2,6) 

58  (0.4) 

70  (0.3) 

95  (0.3) 

6  =  (3,2) 

15  (0.6) 

73  (0.4) 

7  =  (3,7) 

65  (0.5) 

74  (0.4) 

75  (0.1) 

8  =  (3,8) 

59  (0.6) 

72  (0.3) 

98  (0.1) 

9  =  (4,3) 

21  (0.3) 

32  (0.3) 

85  (0.2) 

98  (0.2) 

10  =  (4,9) 

89  (0.7) 

96  (0.3) 

11  =  (5,7) 

32  (0.6) 

48  (0.2) 

67  (0.2) 

12  =  (5,10) 

63  (0.5) 

99  (0.5) 

13  =  (6,3) 

66  (0.8) 

85  (0.1) 

98  (0.1) 

14  =  (6,5) 

6  (0.4) 

15  (0.3) 

39  (0.2) 

58  (0.1) 

15  =  (6,7) 

2  (0.6) 

48  (0.4) 

16  =  (7,8) 

16  (0.3) 

18  (0.3) 

40  (0.2) 

52  (0.2) 

17  =  (7,10) 

16  (0.5) 

34  (0.4) 

71  (0.1) 

18  =  (8,4) 

90  (0.5) 

96  (0.5) 

19  =  (8,9) 

17  (0.4) 

49  (0.4) 

53  (0.1) 

65  (0.1) 

20  =  (9,7) 

3  (0.5) 

30  (0.4) 

50  (0.1) 

21  =  (9,10) 

6  (0.5) 

12  (0.3) 

54  (0.1) 

66  (0.1) 

better  corresponding  GSD  routines.  MCDIST  suffered  badly  on  this  example.  After 
decomposing  10,000,000  sets,  the  probability  of  the  remaining  undetermined  sets 
exceeded  0.80.  This  is  due  to  the  tremendous  size  of  the  state  space  (as  mentioned 
above)  and  to  the  nature  of  the  input  distributions. 

Example  4.3  Taking  the  input  distribution  for  Example  4.2,  we  make  two  changes 
to  make  it  more  reasonable.  First,  we  bring  the  possible  values  for  each  arc  cost 
and  arc  capacity  random  variable  closer  together.  Next,  we  adjust  the  probabilities 
so  that  in  each  case  the  status  quo  is  more  heavily  weighted  than  before.  The 
cardinality  of  the  state  space  remains  the  same  as  that  of  Example  4.2.  The  resulting 
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Table  4.2:  Arc  Capacity  Distributions  for  Example  4.2. 


Arc 

Capacity  (Probability) 

mmm\  jr  n| 

15  (0.7) 

Jj  ■ 

7  (0.2) 

10  (0.2) 

15  (0.6) 

r  ■ 

8  (0.5) 

13  (0.5) 

Hola  H  m 

6  (0.2) 

10  (0.3) 

17  (0.5) 

5  -  (2,6) 

4  (0.1) 

7  (0.2) 

8  (0.3) 

10  (0.4) 

6  =  (3,2) 

7(0.1) 

9  (0.2) 

15  (0.7) 

7  =  (3,7) 

5(0.1) 

8(0.1) 

10  (0.3) 

16  (0.5) 

8  =  (3,8) 

8(0.1) 

12  (0.9) 

9  =  (4,3) 

9(0.1) 

13  (0.2) 

18  (0.7) 

10  =  (4,9) 

4  (0.2) 

7  (0.2) 

11  (0.2) 

113  (0.4) 

11  =  (5,7) 

14  (0.5) 

12  =  (5,10) 

6(0.1) 

15  (0.2) 

17  (0.7) 

13  =  (6,3) 

4(0.1) 

7  (0.2) 

12  (0.2) 

20  (0.5) 

14  =  (6,5) 

5  (0.3) 

12  (0.7) 

15  =  (6,7) 

5  (0.2) 

7  (0.3) 

8  (0.5) 

16  =  (7,8) 

5(1.0) 

17  =  (7,10) 

6  (0.2) 

14  (0.8) 

18  =  (8,4) 

3(0.1) 

5(0.1) 

9  (0.2) 

13  (0.6) 

19  =  (8,9) 

8  (0.2) 

14  (0.6) 

20  =  (9,7) 

4(0.1) 

7  (0.2) 

9  (0.3) 

12  (0.4) 

21  =  (9,10) 

5  (0.3) 

12  (0.7) 
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Table  4.3:  Computing  F(1900)  =  0.3037859483  for  Example  4.2. 


Routine 

Full 

Decompo 

Sets 

sition 

Time 

10,000  Sets  D< 
Bounds 

_ 1 

icompof 

Time 

sed 

\U\ 

P{U) 

MCDA-I 
-Feasible  1 

78,488 

310.25 

0.29394898 

0.31171238 

40.26 

4181 

.0178 

-Feasible  2 

225,270 

699.59 

0.28139953 

0.33257054 

25.40 

2500 

.0512 

-Infeasible  1 

>  500,000 

0.00000000 

0.64673483 

45.44 

1925 

.6467 

-Infeasible  2 

>  500,000 

0.00000000 

0.64673483 

22.63 

1924 

.6467 

-Hybrid 

>  500,000 

0.21396513 

0.47207714 

33.32 

2525 

.2581 

MCDA-II 
-Feasible  1 

>  500,000 

0.09918489 

0.51638467 

28.19 

4463 

.4172 

-Feasible  2 

>  500,000 

0.04420176 

0.87545531 

19.28 

4513 

.8313 

-Infeasible  1 

>  500,000 

0.13529266 

0.39088202 

26.82 

3744 

.2556 

-Infeasible  2 

>  500,000 

0.04796162 

0.40547247 

15.94 

3727 

.3566 

-Hybrid 

>  500,000 

0.12020747 

0.54014936 

16.91 

3786 

.4200 

Table  4.4:  Computing  F( 3500)  =  0.9999996197  for  Example  4.2. 


10,000  Sets  Decomposed 

Bounds 

Time 

\U\ 

P{U) 

0.99999805 

0.99999989 

92.07 

3411 

.000002 

0.99999798 

0.99999990 

77.16 

3413 

.000002 

0.99998046 

1.00000000 

86.03 

4369 

.000020 

0.99768438 

1.00000000 

97.59 

3696 

.002316 

0.99940423 

1.00000000 

71.06 

2762 

.000596 

0.99998638 

1.00000000 

59.16 

2533 

.000014 

0.99998638 

1.00000000 

44.53 

2530 

.000014 

0.99999060 

0.99999993 

57.81 

4725 

.000009 

0.99763126 

0.99999999 

39.88 

3265 

.002369 

0.99965996 

0.99999999 

39.34 

3381 

.000340 

100 


distribution  is  listed  in  Tables  4.5  and  4.6.  The  possible  optimal  flow  costs  now 
range  from  1591  to  2664.  We  evaluate  F(1800)  in  Table  4.7.  Similar  results  for  the 
evaluation  of  F(2000)  are  found  in  Table  4.8.  All  ten  routines  required  more  than 
500,000  iterations  for  the  exact  evaluation  of  F(2000),  too  many  to  be  of  practical 
use.  However,  the  bounds  very  quickly  achieved  a  useable  width.  Despite  the  change 
in  the  data  file,  MCDIST  still  labored.  After  1  million  iterations,  the  probability  of 
the  remaining  undetrmined  sets  was  0.35.  After  10  million,  0.33  still  remained.  A 
partial  listing  of  the  result  of  10  million  iterations  is  given  in  Table  4.9.  Notice  that 
bounds  for  small  values  of  d  are  relatively  tight  as  compared  with  those  for  larger 
d.  This  is  due  to  the  fact  that  procedure  MCDIST  begins  at  the  lowest  possible 
MCF  cost  and  moves  up.  Bounds  on  the  mean  optimal  flow  cost,  calculated  from 
the  bounds  on  this  distribution,  are  1735.4  to  1842.3. 

Example  4.4  Here  the  cost  and  capacity  of  each  arc  form  a  dependent  pair,  as 
shown  in  Table  4.10.  The  cardinality  of  the  state  space  is  about  8  billion.  The 
range  of  optimal  flow  costs  is  from  1591  to  3607.  In  Table  4.11  are  results  for 
F(1900).  The  hybrid  routines  clearly  required  more  computational  effort  than  the 
best  GSD  routines.  Table  4.12  contains  results  for  F(3000)  while  Table  4.13  contains 
a  partial  listing  of  the  cdf  F(d).  Full  decomposition  considered  1,128,522  sets.  The 
mean  optimal  flow  cost  computed  from  this  distribution  is  1979.05,  with  standard 
deviation  251.9. 

Example  4.5  Just  as  we  turned  Example  4.2  into  Example  4.3,  here  we  have  taken 
the  data  from  Example  4.4  and  created  the  present  example.  Again,  the  spread 
of  arc  costs  and  capacities  was  reduced  and  the  probabilities  were  shifted  so  that 
significantly  more  weight  is  on  the  status  quo.  The  resulting  distribution,  given  in 
Table  4.14,  gives  a  range  of  possible  optimal  flow  costs  from  1591  to  2586.  In  Table 
4.15  we  summarize  results  for  evaluating  F(1650).  Results  for  F(2000)  are  in  Table 
4.16.  Again,  the  hybrid  routines  took  siginiflcantly  longer  to  converge  than  the  other 


Table  4.5:  Arc  Cost  Distributions  for  Example  4.3. 


Arc 

Cost  (Probability) 

1  =  (1,2) 

i 

84  (0.1) 

2  =  (1,3) 

1 

IS! 

42  (0.05) 

3  =  (1,4) 

1 

48  (0.1) 

51  (0.1) 

4  =  (2,5) 

26  (0.7) 

31  (0.1) 

45  (0.1) 

50  (0.1) 

5  =  (2,6) 

58  (0.8) 

70  (0.1) 

85  (0.1) 

6  =  (3,2) 

15  (0.9) 

23  (0.1) 

7  =  (3,7) 

65  (0.7) 

74  (0.2) 

75  (0.1) 

8  =  (3,8) 

59  (0.8) 

65  (0.1) 

78  (0.1) 

9  =  (4,3) 

21  (0.7) 

32  (0.1) 

45  (0.1) 

50  (0.1) 

10  =  (4,9) 

89  (0.8) 

96  (0.2) 

11  =  (5,7) 

32  (0.8) 

48  (0.1) 

57  (0.1) 

12  =  (5,10) 

63  (0.8) 

79  (0.2) 

13  =  (6,3) 

66  (0.8) 

75  (0.1) 

88  (0.1) 

14  =  (6,5) 

6  (0.6) 

15  (0.2) 

19  (0.1) 

28  (0.1) 

15  =  (6,7) 

2  (0.8) 

8  (0.2) 

16  =  (7,8) 

16  (0.7) 

18  (0.1) 

22  (0.1) 

35  (0.1) 

17  =  (7,10) 

16  (0.8) 

24  (0.1) 

31  (0.1) 

18  =  (8,4) 

90  (0.9) 

96  (0.1) 

19  =  (8,9) 

17  (0.7) 

29  (0.1) 

33  (0.1) 

45  (0.1) 

20  =  (9,7) 

3  (0.7) 

10  (0.2) 

15  (0.1) 

21  =  (9,10) 

6  (0.6) 

12  (0.2) 

24  (0.1) 

26  (0.1) 

Table  4.6;  Arc  Capacity  Distributions  for  Example  4.3. 


Arc 

Capacity  (Probability) 

■Ijlftfll 

iaroEiJi 

l|||9  B  m 

msm 

wml 

15  (0.7) 

vis  IDS 

8  (0.2) 

13  (0.8) 

BBS 

10  (0.1) 

mvfli 

5  =  (2,6) 

S|CT  1 

7  (0.1) 

Kimi 

in  (0.6) 

6  =  (3,2) 

9  (0.2) 

15  (0.7) 

7  ■-=  (3,7) 

5  (0.1) 

8  (0.1) 

10  (0.1) 

16  (0.7) 

8  =  (3,8) 

8  (0.1) 

12  (0.9) 

9  =  (4,3) 

9  (0.1) 

13  (0.2) 

18  (0.7) 

10  =  (4,9) 

4  (0.1) 

7  (0.1) 

11  (0.2) 

13  (0.6) 

11  =  (5,7) 

8  (0.1) 

14  (0.9) 

12  =  (5,10) 

6  (0.1) 

15  (0.2) 

17  (0.7) 

13  =  (6,3) 

4  (0.1) 

7  (0.1) 

12  (0.1) 

18  (0.7) 

14  =  (6,5) 

5(0.1) 

12  (0.9) 

15  =  (6,7) 

5  (0.1) 

7  (0.1) 

8  (0.8) 

16  =  (7,8) 

5  (1.0) 

17  =  (7,10) 

6  (0.2) 

14  (0.8) 

18  =  (8,4) 

3  (0.1) 

5  (0.1) 

9  (0.2) 

13  (0.6) 

19  =  (8,9) 

8  (0.1) 

10  (0.2) 

14  (0.7) 

20  =  (9,7) 

4  (0.1) 

7  (0.1) 

9  (0.2) 

12  (0.6) 

21  =  (9,10) 

5  (0.1) 

12  (0.9) 
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Table  4.7:  Computing  F(1800)  =  0.6675061884  for  Example  4.3. 


1 

Routine 

Full 

Decompo 

Sets 

sition 

Time 

10,000  Sets  D< 
Bounds  j 

ccompcx 

Time 

sed 

P\ 

P{U) 

MCDA-I 
-Feasible  l 

21,148 

76.78 

0.66744115 

0.66752595 

40.66 

3840 

.00009 

-Feasible  2 

59,252 

130.68 

0.66412553 

0.66934024 

25.91 

2362 

.00522 

-Infeasible  1 

>  500,000 

0.01365559 

0.90549725 

64.88 

2259 

.89184 

-Infeasible  2 

>  500,000 

0.010605*0 

0.90549840 

51.06 

2283 

.89489 

-Hybrid 

0.62086991 

0.69095869 

33.41 

2564 

.07009 

MCDA-II 
-Feasible  l 

>  500,000 

0.51748501 

0.69914193 

25.06 

4258 

.18166 

-Feasible  2 

>  500,000 

0.46837277 

0.80181865 

19.38 

4000 

.33340 

-Infeasible  1 

335,712 

783.19 

0.59588386 

0.68162250 

26.38 

3350 

.08574 

-infeasible  2 

>  500,000 

0.43849046 

0.69493580 

17.28 

3367 

.25640 

-Hybrid 

>  500,000 

0.52844341 

0.72242710 

19.09 

3251 

.19398 

Table  4.8:  Bounding  F(2000)  for  Example  4.3. 


Routine 

10,000  Sets  D 
Bounds 

ecompo 

Time 

sed 

w 

P{U) 

MCDA-I 
-Feasible  1 

0.91922537 

0.97045163 

40.16 

4751 

.05123 

-Feasible  2 

0.91259465 

0.97513833 

27.85 

4219 

.06254 

-Infeasible  1 

0.05918410 

0.99954251 

31.44 

2678 

.94036 

-Infeasible  2 

0.00790706 

0.99955266 

24.57 

2794 

.99165 

-Hybrid 

0.81577920 

0.99071067 

33.47 

3294 

.17493 

MCDA-II 
-Feasible  1 

0.79366799 

0.97677949 

26.34 

4643 

.18311 

-Feasible  2 

0.78470934 

0.99022890 

18.59 

4591 

.20552 

-Infeasible  1 

0.51684255 

0.9838851 

27.62 

4715 

.46704 

-Infeasible  2 

0.29276291 

0.98467936 

18.41 

4379 

.69192 

-Hybrid 

0.77839503 

0.98117291 

19.35 

3760 

.20280 
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Tabic  4.9:  Bounding  the  cdf  F(d)  for  Example  4.3. 


V 

a 

W) 

F.(d) 

d 

1 

1 

mmm 

1600 

0.1390231 

0.1390231 

2150 

0.8421683 

0.9998394 

1650 

0.3164801 

0.3274079 

2200 

0.8623073 

0.9999785 

1700 

0.3775094 

0.4164037 

2250 

0.8903911 

0.9999958 

1750 

0.4869543 

0.5755966 

2300 

0.9221442 

0.9999993 

1800 

0.5640803 

0.7176557 

2350 

0.9459363 

0.9999999 

1850 

0.6247240 

0.8149941 

2400 

0.9685088 

0.9999999 

1900 

0.6932091 

0.9173447 

2450 

0.9813164 

0.9999999 

1950 

0.7380265 

0.9631111 

2500 

0.9906127 

1.0000000 

2000 

0.7690676 

0.9893741 

2550 

0.9965654 

1.0000000 

2050 

0.7918757 

0.9961165 

2600 

0.9983520 

1.0000000 

2100 

0.8171786 

0.9931391 

2650 

0.9998800 

1.0000000 

Table  4.10:  Input  Probability  Distribution  for  Example  4.4. 


1)  (94,10)  (0.2) 

l)  (42,7)  (0.1) 

1)  (61,8)  (0.2) 

1)  (55,8)  (0.2) 

1)  (95,7)  (0.3) 

l\ 

(88,6)  (0.1) 

V 

\)  (75.8)  (0.1) 

0  (98,8)  (0.1) 

1)  (85,12)  (0.2) 

(98,10)  (0.2) 

!)  (67,8)  (0.2) 

)  (98,7)  (0.1) 
t)  (39,8)  (0.2) 

(58,5)  (0.1) 

1)  (40,8)  (0.2) 

1)  (71,6)  (0.1) 

(52,5)  (0.2) 

1)  (53,10)  (0.1) 

1)  (50,7)  (0.1) 

1)  (54,8)  (0.1) 

(65,8)  (0.1) 

(66,5)  (0.1) 
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Table  4.11:  Computing  F(1900)  =  0.36696615  for  Example  4.4. 


Routine 

Full 

Decompo 

Intervals 

sition 

Time 

500  Sets  De 
Bounds 

compost 

Time 

ed 

\U\ 

P{U) 

MCDA-I 
-Feasible  1 

559 

3.81 

0.36692783 

0.36699223 

3.38 

32 

.000065 

-Feasible  2 

1497 

6.16 

0.36368208 

0.37200554 

2.43 

135 

.008323 

-Infeasible  1 

3179 

22.09 

0.34150064 

0.37240647 

6.22 

248 

.030906 

-Infeasible  2 

3179 

22.20 

0.34150064 

0.37240647 

6.34 

248 

.030906 

-Hybrid 

0.23197761 

0.39739952 

13.32 

1662 

.165422 

MCDA-II 
-Feasible  1 

>  100,000 

0.29763540 

0.41543235 

6.84 

2274 

.1178 

-Feasible  2 

>  100,000 

0.24929387 

0.65997558 

5.46 

.4107 

-Infeasible  1 

21,469 

124.15 

0.30655235 

0.41347403 

5.56 

1114 

.10695 

-Infeasible  2 

>  100, Ooo 

0.24629346 

0.42723007 

4.18 

1413 

.18094 

-Hybrid 

59,857 

224.80 

C.32552317 

0.40912128 

6.06 

1786 

.083598 
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Table  4.12:  Computing  F(3000)  =  0.99908320  for  Example  4.4. 


Routine 

Full 

Decompo 

Intervals 

sition 

Time 

500  Sets  De 
Bounds 

compost 

Time 

ed 

P\ 

P{U) 

MCDA-I 
-Feasible  1 

2908 

24.06 

0.99875743 

0.99930461 

4.84 

445 

.000547 

-Feasible  2 

6127 

30.78 

0.99869394 

0.99940407 

4.79 

412 

.000710 

-Infeasible  1 

8455 

42.03 

0.97500671 

0.99931205 

5.50 

405 

.024305 

-Infeasible  2 

8455 

42.72 

0.97500671 

0.99931205 

5.45 

405 

.024305 

-Hybrid 

71,439 

299.40 

0.95753799 

0.99949807 

8.20 

709 

.041960 

MCDA-II 
-Feasible  1 

>  100,000 

0.99793109 

0.99958227 

5.41 

850 

.001651 

-Feasible  2 

>  100,000 

0.99789501 

0.99970247 

4.59 

827 

.001808 

-Infeasible  1 

19,711 

108.85 

0.99794290 

0.99950216 

6.44 

.001559 

-Infeasible  2 

>  100,000 

0.83254224 

0.99997934 

5.34 

3351 

.167440 

-Hybrid 

55,659 

215.15 

0.99809057 

0.99945495 

5.87 

.001364 

Table  4.13:  Partial  listing  of  the  cdf  F(d)  for  Example  4, 


d 

F(d ) 

d 

m 

1600 

0.0240000 

2150 

0.7852381 

1650 

0.1285435 

2200 

0.8358933 

1700 

0.1665479 

2250 

0.8707846 

1750 

0.2019579 

2300 

0.9006702 

1800 

0.2400117 

2350 

0.9241992 

1850 

0.2834110 

2400 

0.9469183 

1900 

0.3669662 

2500 

0.9722436 

1950 

0.4770354 

2600 

0.9821508 

2000 

0.5725981 

2700 

0.9887615 

2050 

0.6538413 

3000 

0.9990832 

2100 

0.7296389 

3500 

0.9999979 

Table  4.14:  Input  Probability  Distribution  for  Example  4.5. 


Arc 

(Capacity, Cost' 

(Probability) 

1  =  (1,2) 

(70,15)  (0.7) 

(73,12)  (0.2) 

(84,10)  (0.1) 

2  --  (1,3) 

(25,15)  (0.9) 

(35,10)  (0.05) 

(42,7)  (0.05) 

3  =  (1,4) 

(42,13)  (0.8) 

(48,10)  (0.1) 

(51,8)  (0.1) 

4  =  (2,5) 

(26,12)  (0.6) 

(31,10)  (0.2) 

(45,8)  (0.1) 

(50,6)  (0.1) 

5  =  (2,6) 

(58,10)  (0.8) 

(70,8)  (0.1) 

(85,7)  (0.1) 

6  =  (3,2) 

(15,15)  (0.9) 

(23,9)  (0.1) 

7  =  (3,7) 

(65,16)  (0.7) 

(74,10)  (0.2) 

(75,8)  (0.1) 

8  =  (3,8) 

(59,12)  (0.8) 

(65,10)  (0.1) 

(78,8)  (0.1) 

9  =  (4,3) 

(21,18)  (0.6) 

(32,13)  (0.2) 

(45,12)  (0.1) 

(50,10)  (0.1) 

10  =  (4,9) 

(89,11)  (0.8) 

(96,7)  (0.2) 

11  =  (5,7) 

(32,14)  (0.7) 

(48,12)  (0.2) 

(57,8)  (0.1) 

12  =  (5,10) 

(63,17)  (0.8) 

(79,15)  (0.2) 

13  =  (6,3) 

(G6,18)  (0.8) 

(75,12)  (0.1) 

(88,7)  (0.1) 

14  =  (6,5) 

(6,12)  (0.7) 

(15,10)  (0.1) 

(19,8)  (0.1) 

(28,5)  (0.1) 

15  =  (6,7) 

(2,8)  (0.9) 

(8,7)  (0.1) 

16  =  (7,81 

(16,12)  (0.6) 

(18,11)  (0.2) 

(22,8)  (0.1) 

(35,5)  (0.1) 

17  =  (7,10) 

(16,14)  (0.8) 

(24,10)  (0.1) 

(31,6)  (0.1) 

18  =  (8,4) 

(90,13)  (0.8) 

(96,9)  (0.2) 

19  =  (8,9) 

(17,14)  (0.7) 

(29,12)  (0.1) 

(33,10)  (0.1) 

(45,8)  (0.1) 

20  =  (9,7) 

(3,12)  (0.8) 

(10,9)  (0.1) 

(15,7)  (0.1) 

21  =  (9,10) 

(6,12)  (0.6) 

(12,10)  (0.2) 

(24,8)  (0.1) 

(26,5)  (0.1) 

methods.  Finally,  a  partial  listing  of  the  distribution  F{d)  is  given  in  Table  4.17. 
Full  decomposition  considered  448,729  intervals.  The  mean  optimal  flow  cost,  as 
computed  from  this  distribution,  is  1689.7094  with  standard  deviation  135.01. 


4.5  Conclusions 

We  conclude  this  chapter  by  making  some  observations  about  the  use  of  parti¬ 
tioning  routines  to  solve  the  stochastic  minimum  cost  flow  problem.  Several  points 
deserve  mentioning  at  this  time. 
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Table  4.15:  Computing  F(1650)  =  0.62234033  for  Example  4.5. 


Full 

500  Sets  Decomposed 

Routine 

Intervals 

Time 

Bounds 

Time 

\U\ 

P{U} 

MCDA-I 

-4# 

-Feasible  1 

15 

1 

-Feasible  2 

30 

| m 

-Infeasible  1 

30 

ia 

-Infeasible  2 

30 

m 

-Hybrid 

3,514 

14.56 

0.61302130 

0.62792051 

5.13 

156 

.014899 

MCDA-II 

-Feasible  1 

775 

3.53 

0.62234028 

0.62234033 

2.60 

5.2E-8 

-Feasible  2 

>  100,000 

0.53768333 

0.76617651 

5.12 

.22849 

-Infeasible  1 

28 

0.53 

-Infeasible  2 

6011 

19.09 

0.62233192 

0.62234056 

2.06 

1046 

8.6E-6 

-Hybrid 

1094 

5.13 

0.62232657 

0.62234094 

3.16 

47 

1.4E-5 

lil 
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Table  4.16:  Computing  F(2000)  =  0.95451624  for  Example  4.5. 


Routine 


MCDA-I 
-Feasible  1 
-Feasible  2 
-Infeasible  1 
-Infeasible  2 
-Hybrid 


MCDA-II 
-Feasible  1 
-Feasible  2 
-Infeasible  1 
-Infeasible  2 
-Hybrid 


Full 

Decomposition 
Intervals  |  Time 


10.68 

16.56 

21.31 

20.63 


>  100,000 


>  100,000 
>  100,000 
19,346 
>  100,000 
61,931 


116.28 


264.53 


500  Sets  Decomposed 
Bounds  |  Time  |  \U\ 

pm 

0.95414827 

0.95469526 

4.56 

321 

.000547 

0.95362469 

0.95519911 

3.41 

282 

.001574 

0.95074918 

0.95516215 

5.84 

190 

.004413 

0.95074918 

0.95516215 

5.44 

190 

.004413 

0.85767863 

0.96448481 

8.32 

1037 

.106806 

0.94053428 

0.95723072 

5.54 

1491 

.016696 

0.93344046 

0.96810572 

4.87 

2283 

.034665 

0.94785747 

0.95657733 

5.97 

1175 

.008720 

0.81657072 

0.96065541 

4.28 

1977 

.144085 

0.94789470 

0.95614263 

6.35 

1260 

.008248 

Table  4.17:  Partial  listing  of  cdf  F(d)  for  Example  4.5. 


d 

F(d) 

d 

F(d) 

1600 

0.2469600 

1850 

0.8618947 

1625 

0.4921901 

1900 

0.9090322 

1650 

0.6223403 

2000 

0.9545162 

1675 

0.6355670 

2100 

0.9916873 

1700 

0.6593318 

2200 

0.9982173 

1725 

0.7007062 

2300 

0.9997193 

1750 

0.7725537 

2400 

0.9999767 

1775 

0.8127293 

2500 

0.9999996 

1800 

0.8348261 

2600 

1.0000000 

1825  0.8431940 


•  It  is  clear  that,  compared  to  the  success  of  the  MST  routines  in  Chapter  3, 
the  results  of  the  present  chapter  do  not  stand  out  as  being  phenomenal. 

•  Since  each  arc  now  has  two  associated  random  variables  the  size  of  the  state 
space  is  effectively  squared.  Taking  this  into  consideration,  most  of  the  decom¬ 
position  results  here  perform  well  in  terms  of  only  evaluating  a  small  fraction 
of  the  points  in  the  state  space. 

•  It  seems  clear  that  our  conjectures  in  Section  3.5  remain  valid  for  MCF  prob¬ 
lems.  Since  these  problems  do  not  have  a  matroid  structure,  we  did  not  expect 
the  kind  of  performance  we  saw  in  Chapter  3.  In  particular,  the  lack  of  an 
easy  way  to  determine  a  new  minimum  cost  flow  when  arc  costs  or  capacities 
changed  really  hurt.  We  were  left  either  using  a  cutoff  state  that  was  not  deep 
enough  or  reoptimizing  over  and  over  to  try  to  improve  such  a  point. 

•  In  Chapter  3,  the  simpler  routine  STDA-II  (which  did  not  derive  infeasible 
sets  by  finding  limiting  indices)  was  more  effective  than  the  more  complete, 
sophisticated  STDA-I.  Here  we  have  the  opposite  case.  The  inability  to  derive 
easy  cut-off  states  necessitates  the  derivation  of  infeasible  sets. 

•  Here  it  is  reasonable  only  to  compare  the  hybrid  routines  to  the  GSD  routines 
that  are  based  on  the  same  cutoff  state  derivations  (i.e.,  MCDA-I  or  MCDA- 
II).  In  all  cases,  as  predicted  in  Section  2.5,  these  did  not  perform  as  well  as 
the  better  GSD  routines. 

In  summary,  the  algorithms  given  here,  especially  in  light  of  the  above  comments, 
have  performed  very  well.  The  astronomically  large  state  spaces  and  lack  of  special 
matroid  structure  certainly  hindered  their  effectiveness  somewhat.  But  considering 
the  impossibly  hard  task  of  evaluating  these  measures  in  the  absence  of  such  routines 
(by  enumeration  or  crude  Monte  Carlo  simulation),  the  algorithms  of  this  chapter 
must  be  called  successful. 


Finally,  there  are  several  important  network  optimization  problems  that  can  be 
cast  as  minimum  cost  network  flow  problems.  Furthermore,  independence  among 
components  in  these  problems  translates  into  independence  in  the  present  setting  so 
that  application  of  our  methods  makes  sense.  However,  we  point  out  that  we  would 
expect  for  GSD-type  algorithms  tailored  to  exploit  the  special  structure  of  one  of 
these  problems  to  out-perform  the  more  general  algorithms  given  here  when  applied 
to  the  same  problem.  Indeed,  this  has  been  seen  to  be  the  case  [2]. 
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CHAPTER  5 

Conclusions  and  Directions  for 
Further  Study 


» 
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In  the  preceding  chapters  we  have  attempted  to  accomplish  several  goals.  First,  } 

we  have  tried  to  improve  the  state  of  knowledge  concerning  solution  methodologies 
that  use  iterative  partitioning  of  state  spaces.  Secondly,  we  have  endeavored  to 
adapt  these  procedures  to  create  efficient  means  for  studying  minimum  spanning 
trees  in  stochastic  networks.  Finally,  we  have  attempted  to  likewise  address  certain 
probabilistic  minimum  cost  network  flow  problems.  We  discuss  here  the  extent  to 
which  we  have  succeeded  in  these  aims,  and  propose  directions  for  further  study. 

In  Chapter  2,  our  task  was  to  investigate  a  computational  approach  that  has  been  | 

used  with  some  success  in  addressing  other  problems  involving  stochastic  network 
systems.  This  approach,  best  described  as  a  variant  of  factoring,  was  shown  to 
be  applicable  to  a  rather  broad  class  of  stochastic  systems.  Our  description  of  the 
General  State  Space  Decomposition  (GSD)  algorithm  was  intended  in  part  to  make  a 
contribution  by  being  complete  and  flexible  enough  to  make  this  powerful  technique 
easily  adaptable  to  all  problems  in  this  class  of  systems. 

In  addition,  our  description  was  designed  to  facilitate  some  theoretical  investi-  ) 

gations.  Specifically,  we  posed  some  questions  concerning  the  partitioning  step  used 
in  GSD  in  hopes  that  their  answers  might  shed  some  light  on  certain  performance 
issues.  We  essentially  sought  to  find  the  best  approach  to  removing  feasible  and 
infeasible  states  from  an  undetermined  interval.  Our  analysis  was  based  on  the 
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intuitive  notion  that  performance  would  suffer  whenever  an  interval  was  split  into 
many  small  pieces  instead  of  fewer,  larger  pieces.  Our  conclusion,  that  any  depar¬ 
ture  from  GSD  would  result  in  an  unnecessary  proliferation  of  filed  undetermined 
sets,  led  us  to  predict  that,  in  general,  the  better  GSD  routines  should  out-perform 
arbitrary  decomposition  schemes  applied  to  the  same  problem.  This  conjecture  was 
supported  by  virtually  every  application  studied  in  Chapters  3  and  4. 

Nevertheless,  there  are  some  performance  issues  that  remain  to  be  investigated. 
These  difficult  issues  deal  with  finding  the  best  GSD  scheme  for  a  given  problem. 
Why,  for  example,  do  certain  cutoff  state  derivations  perform  well  on  some  instances 
and  poorly  on  others?  Given  a  number  of  alternative  cutoff  state  derivations,  is  there 
a  way  to  determine  which  derivation  will  perform  best  for  a  particular  problem 
instance?  These  questions  are  difficult  because  of  the  complex  interplay  of  different 
factors  that  determine  a  GSD  algorithm’s  ultimate  effectiveness  for  a  given  problem, 
but  are  nonetheless  worthy  of  further  exploration. 

The  subject  of  Chapter  3  was  the  application  of  GSD  to  problems  related  to  min¬ 
imum  spanning  trees  in  stochastic  networks.  Both  problems  addressed  there  were 
shown  to  be  computationally  hard,  a  heretofore  unproved  result  for  the  criticality 
index  problem.  Despite  this,  we  found  that  these  algorithms  were  very  effective. 
We  hypothesized  that  this  was  due  in  part  to  the  matroid  structure  of  the  MST 
problem.  The  fact  that  the  algorithms  in  Chapter  4  (applied  to  a  problem  with  no 
matroid  structure)  performed  worse  seemed  to  support  this  conjecture. 

The  routine  used  to  compute  criticality  indices  stood  out,  as  being  particularly 
efficient.  Motivated  by  this  faOt,  we  proposed  in  Section  3.6  a  heuristic  algorithm 
for  determining  the  membership  of  the  most  probable  minimum  spanning  tree.  We 
believe  this  to  be  worthy  of  study  and  intend  to  investigate  this  application  in 
the  future.'  In  addition,  from  a  more  theoretical  standpoint,  it  might  be  fruitful 
to  look  for  underlying  structural  reasons  for  the  success  of  GSD  applied  to  this 
problem.  Such  information  might  help  us  to  predict  with  improved  confidence  the 
kinds  of  problems  that  will  readily  admit  solution  by  GSD  and  to  better  tailor 
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specific  decomposition  methods  to  those  problems. 

In  Chapter  4  we  applied  partitioning  algorithms  to  stochastic  minimum  cost 
network  flow  problems.  These  problems,  too,  were  shown  to  be  #P-hard.  This 
gave  us  cause  for  concern  since  in  this  setting  each  arc  has  two  associated  random 
variables,  resulting  in  state  spaces  that  are  larger  than  those  for  corresponding  MST 
problems.  Furthermore,  as  mentioned  above,  the  MCF  problem  has  no  matroid 
structure.  Consistent  with  this  expectation,  computational  performance  was  not  as 
as  remarkable  as  that  which  we  saw  in  Chapter  3.  However,  we  were  still  able  to  get 
usable  computational  results  in  most  cases,  and  these  results  continued  to  confirm 
the  results  of  Section  2.5.  Particularly  effective  were  the  algorithms  applied  to  the 
case  of  dependent  arc  costs  and  capacities.  These  applications  were  also  important 
because  they  illustrated  a  new  flexibility  of  GSD  and  represented,  as  far  as  we  know, 
the  first  time  vector-valued  variables  were  used  in  a  decomposition  approach  of  this 
kind. 

However,  for  large  problems  it  would  likely  be  worthwhile  to  write  variance- 
reducing  Monte  Carlo  estimation  routines  of  the  form  given  in  Section  2.7.  This 
investment  in  additional  effort  would  allow  us  to  get  a  good  set  of  initial  bounds 
for  PF  via  GSD,  and  then  use  the  remaining  undetermined  intervals  to  efficiently 
produce  an  unbiased  estimate  of  PF. 

This  thesis  has  made  a  number  of  contributions.  First,  it  has  shown  several 
important  stochastic  network  problems  to  be  computationally  hard,  and  then  has 
presented  efficient  means  for  solving  these  problems  exactly.  In  addition,  it  has 
advanced  the  understanding  of  state  space  partitioning  methods.  In  doing  so,  it 
has  made  these  methods  more  accessible  and  has  drawn  some  strong  conclusions 
about  the  soundness  of  using  GSD-type  routines  instead  of  arbitrary  decomposition 
schemes.  Finally,  it  haw  proposed  some  areas  in  which  further  gains  might  be  made 
with  regards  to  this  powerful  computational  technique. 


I 


» 


I 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


4* 


118 


ft 


I 

Bibliography 

ft 

[1]  Alexopoulos,  C.  (1991).  Computing  criticality  indices  of  arcs  and  the  mean 

maximum  flow  value  in  networks  with  discrete  random  capacities,  Technical  ) 

Report,  School  of  Industrial  and  Systems  Engineering,  Georgia  Institute  of 
Technology,  submitted  for  publication. 

[2]  Alexopoulos,  C.  (1992).  State  space  partitioning  methods  for  stochastic  shortest  ft 

path  problems,  Technical  Report,  School  of  Industrial  and  Systems  Engineering, 

Georgia  Institute  of  Technology. 

[3]  Alexopoulos,  C.  and  G.S.  Fishman  (1991).  Characterizing  stochastic  flow  net-  ft 

works  using  the  Monte  Carlo  method,  Networks,  21,  775-798. 

[4]  Alexopoulos,  C.  and  G.S.  Fishman  (1991).  Sensitivity  analysis  in  stochastic 

flow  networks  using  the  Monte  Carlo  Method,  Networks,  forthcoming.  ) 

[5]  Bailey,  M.P.  (1990).  Minimization  on  stochastic  matroids,  NPS-55-90-01,  Naval 
Postgraduate  School,  Monterrey,  California. 

[6]  Ball,  M.O.  (1987).  Computational  complexity  of  network  reliability  analysis: 

An  overview,  IEEE  Transactions  on  Reliability  35,  230-239. 

[7]  Ball,  M.O.,  Colbourn,  C.J.  and  J.S.  Provan  (1992).  Network  Reliability  (Sub¬ 
mitted  for  publication).  * 

[8]  Barlow,  R.E.  and  Proschan,  F.  (1981).  Statistical  Theory  of  Reliability  and  Life 
Testing ,  Holt,  Reinhart  and  Winston. 

I 


ft 


119 


[9]  Bertsekas,  D.P.  and  P.  Tseng  (1988)  The  RELAX  codes  for  linear  minimum 
cost  network  flow  problems,  Annals  of  Operations  Research ,  7,  125-190. 

[10]  Bertsimas,  D.J.  (1990).  The  probabilistic  minimum  spanning  tree  problem,  Net¬ 
works ,  20,  245-175. 

[11]  Bertsimas,  D.J.  and  G.  van  Rysin  (1990).  An  asymptotic  determination  of 
the  minimum  spanning  tree  and  minimum  matching  constants  in  geometrical 
probability,  Operations  Research  Letter ,  9,  223-231. 

[12]  Camerini,  P.M.,  Galbiati,  G.  and  F.  Maffioli  (1988).  Algorithms  for  finding 
optimum  trees:  description,  use  and  valua  .on,  Annals  of  Operations  Research , 
IS,  265-397. 

[13]  Chen,  H.  and  J.  Zhou  (1991).  Reliability  optimization  in  generalized  stochastic- 
flow  networks,  IEEE  Transactions  on  Reliability ,  40,  92-97. 

[14]  Corea,  G.A.  and  V.G.  Kulkarni  (1990).  Minimum  cost  routing  on  stochastic 
networks,  Operations  Research ,  38,  527-536. 

[15]  Doulliez,  P.  and  E.  Jamoulle  (1972).  Transportation  networks  with  random  arc 
capacities,  R.A.I.R.O. ,  3,  45-60. 

[16]  Fishman,  G.S.  and  T.  Shaw  (1989).  Evaluating  reliability  of  stochastic  flow 
networks,  Probability  in  the  Engineering  and  Informational  Sciences ,  3,  493- 
509. 

[17]  Frieze,  A.M.  (1985).  On  the  value  of  a  random  minimum  spanning  tree  problem, 
Discrete  Applied  Mathematics,  10,  47-56. 

[18]  Gilbert,  E.M.  (1965).  Random  minimal  trees,  J.  SIAM,  13,  376-387. 

[19]  Grimmett,  G.R.  and  D.J. A.  Welsh  (1982).  Flow  in  networks  with  random  ca¬ 
pacities,  Stochastics,  7,  205-229. 


120 


» 

I 

[20]  Hagstrom,  J.N.  (1983).  Using  the  decomposition-tree  of  a  network  in  reliability  ^ 

computation,  IEEE  Transactions  on  Reliability ,  32,  71-78. 

[21]  Hagstrom,  J.N.  and  P.  Kumar  (1984).  Reliability  computation  on  a  probabilistic 

network  with  path  length  criterion,  Technical  Report,  University  of  Illinois,  * 

Chicago. 

[22]  Jain,  A.  and  J.W.  Mamer  (1988)  Approximations  for  the  random  minimal  span¬ 
ning  tree  with  applications  to  network  provisioning,  Operations  Research ,  36,  ^ 

575-584. 

[23]  Kulkarni,  V.G.  (1986).  Minimum  spanning  trees  in  undirected  networks  with 

exponentially  distributed  arc  weights,  Networks,  16,  111-124.  * 

[24]  Lawler,  E.L.  (1976).  Combinatorial  Optimization:  Networks  and  Matroids, 

Holt,  Rinehart  and  Winston. 

ft 

[25]  Le,  K.V.  and  V.O.K.  Li  (1989).  Modeling  and  analysis  of  systems  with  mul¬ 
timode  components  and  dependent  failures,  IEEE  Transactions  on  Reliability , 

29,  68-75. 

ft 

[26]  Leuker,  G.S.  (1981).  Optimisation  Problems  on  graphs  with  independent  edge 
weights,  SIAM  Journal  of  Computing,  10,  338-351. 

[27]  Nagamochi,  H.  and  T.  Ibaraki  (1992).  On  Onaga’s  upper  bound  on  the  mean 
value  of  probabilistic  maximum  flows,  IEEE  Transactions  on  Reliability ,  41, 

225-229. 

[28]  Nemhauser,  G.L.  and  L.A.  Wolsey  (1988).  Integer  and  Combinatorial  Optimiza -  ^ 

tion,  Wiley- Interscience. 

[29]  Papadimitriou,  C.  and  K.  Stieglitz  (1982).  Combinatorial  Optimization:  Algo¬ 
rithms  and  Complexity,  Prentice  Hall. 


» 


121 


> 


[30]  Prekopa,  A.  and  E.  Boros  (1991).  On  the  existence  of  a  feasible  flow  in  a 
stochastic  transshipment  network,  Operations  Research ,  39,  119-129. 

[31]  Steele,  J.M.  (1987)  Growth  rates  on  Frieze’s  £(3)  limit  for  lengths  of  minimal 
spanning  trees,  Discrete  Applied  Mathematics,  18,  99-103. 


ft 


» 


ft 


ft 


ft 


I 


ft 


ft 


122 


VITA 


Captain  Jay  Alan  Jacobson  was  born  to  James  Edmund  and  Diana  Sue 
Jacobson  on  September  15,  1961,  in  Birmingham,  Alabama.  He  attended  Huffman 
High  School.  Upon  graduating  in  1979,  he  entered  the  University  of  Alabama  under 
a  four-year  Air  Force  ROTC  scholarship  to  study  mathematics.  In  1983  he  was 
awarded  the  Bachelor  of  Science  degree  ( magna  cum  laude)  and  was  commissioned 
an  officer  in  the  United  States  Air  Force. 

Air  Force  assignments  have  included  serving  as  Chief  of  Airman  Modeling  at 
the  Air  Force  Military  Personnel  Center  at  Randolph  AFB,  Texas,  and  being  an 
instructor  in  the  Department  of  Mathematical  Sciences  at  the  United  States  Air 
Force  Academy  in  Colorado. 

He  was  awarded  a  Master  of  Science  degree  in  operations  research  from  Stanford 
University  in  1988.  He  began  studies  in  the  School  of  Industrial  and  Systems  Engi¬ 
neering  at  the  Georgia  Institute  of  Technology  in  1989,  and  was  awarded  the  degree 
of  Doctor  of  Philosophy  in  1993. 

He  is  married  to  the  former  Jennifer  Karen  Swann  of  Northport,  Alabama.  They 
have  two  children.  John  Christian  was  born  in  1985  at  Fort  Sam  Houston,  Texas, 
and  William  Alan  was  born  in  1989  at  the  Air  Force  Academy. 


123 


